Skip to content

Commit

Permalink
Merge pull request #167 from JuliaGNI/dev-variational
Browse files Browse the repository at this point in the history
Implement additional variational and Hamilton-Pontryagin integrators
  • Loading branch information
michakraus authored Aug 30, 2023
2 parents 10ed680 + e0132c8 commit 1f0da91
Show file tree
Hide file tree
Showing 28 changed files with 972 additions and 81 deletions.
3 changes: 2 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,8 @@ makedocs(bib,
"DVI" => "integrators/dvi.md",
# "CGVI" => "integrators/cgvi.md",
# "DGVI" => "integrators/dgvi.md",
# "HPG" => "integrators/hpg.md",
"Hamilton-Pontryagin" => "integrators/hpi.md",
# "Hamilton-Pontryagin-Galerkin" => "integrators/hpg.md",
],
"Modules" => [
# "Discontinuities" => "modules/discontinuities.md",
Expand Down
58 changes: 38 additions & 20 deletions docs/src/index.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

# GeometricIntegrators.jl

*Julia library of geometric integrators for differential equations.*
Expand Down Expand Up @@ -47,38 +46,57 @@ While DifferentialEquations.jl provides a feature-rich ecosystem for the solutio
## Manual

```@contents
Pages = ["tutorial/tutorial.md",
"problems.md",
"methods.md",
Pages = [
"tutorial/tutorial.md",
"problems.md",
"methods.md",
]
```


## Modules
## Integrators

```@contents
Pages = ["modules/methods.md",
"modules/integrators.md",
#"modules/discontinuities.md",
"modules/problems.md",
"modules/equations.md",
#"modules/simulations.md",
"modules/solutions.md",
"modules/rungekutta.md",
"modules/rungekutta_partitioned.md",
"modules/spark.md",
Pages = [
"integrators/rk.md",
"integrators/splitting.md",
"integrators/variational.md",
"integrators/vprk.md",
"integrators/spark.md",
"integrators/dvi.md",
"integrators/cgvi.md",
"integrators/dgvi.md",
"integrators/hpi.md",
"integrators/hpg.md",
]
```

## Modules

```@contents
Pages = [
"modules/methods.md",
"modules/integrators.md",
# "modules/discontinuities.md",
"modules/problems.md",
"modules/equations.md",
# "modules/simulations.md",
"modules/solutions.md",
"modules/rungekutta.md",
"modules/rungekutta_partitioned.md",
"modules/spark.md",
]
```

## Developer Documentation

```@contents
Pages = ["developer/integrators.md",
"developer/projections.md",
"developer/code_integration.md",
"developer/custom_integrators.md",
"developer/adaptive_time_stepping.md",
Pages = [
"developer/integrators.md",
"developer/projections.md",
"developer/code_integration.md",
"developer/custom_integrators.md",
"developer/adaptive_time_stepping.md",
]
```

Expand Down
109 changes: 109 additions & 0 deletions docs/src/integrators/hpi.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
# [Hamilton-Pontryagin Integrators](@id hpi)



## Trapezoidal Discrete Lagrangian

We consider a discrete Lagrangian of the form
```math
L_d (q_{n}, q_{n+1}) = \frac{h}{2} \big[ L (q_{n}, v_{n+1/2}) + L (q_{n+1}, v_{n+1/2}) \big] ,
```
where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$.
The discrete Hamilton-Pontryagin action reads
```math
A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] ,
```
where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$.
A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is
```math
\phi_h (q_{n}, q_{n+1}; a_{n,n+1}) = \frac{q_{n+1} - q_{n}}{h} .
```
In order to solve the discrete Euler-Lagrange equations, the user needs to specify the map $\phi_h$, its gradients with respect to $q_{n}$ and $q_{n+1}$, denoted by $D_1 \phi_h$ and $D_2 \phi_h$, respectively, the gradient with respect to the parameters, denoted by $D_a \phi_h$, and an initial set of parameters $a_{0}$.

The equations of motion, that are solved by this integrator, are computed as:
```math
\begin{aligned}
0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2})
+ \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n-1/2}) \\
&+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2}
+ h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\
0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\
v_{n+1/2}
&= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\
p_{n+1/2}
&= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2})
+ \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) .
\end{aligned}
```
Upon defining the momentum
```math
p_{n}
= h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2}
+ \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n-1/2}) ,
```
we can rewrite the equations of motion as
```math
\begin{aligned}
0 &= p_{n}
+ \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2})
+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\
0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\
v_{n+1/2}
&= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\
p_{n+1/2}
&= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2})
+ \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) , \\
p_{n+1}
&= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2}
+ \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n+1}, v_{n+1/2}) .
\end{aligned}
```
Given $(q_{n}, p_{n})$, the first four equations can be solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method, and the last equation provides an explicit update for $p_{n+1}$.


## Midpoint Discrete Lagrangian

Similarly to the previous section, we can consider a discrete Lagrangian of the form
```math
L_d (q_{n}, q_{n+1}) = h \, L \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) ,
```
where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$.
The discrete Hamilton-Pontryagin action reads
```math
A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_{d}^{\alpha,\beta} (q_{n}, q_{n+1}, v_{n+1/2}) + \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] ,
```
where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$.

The equations of motion, that are solved by this integrator, is computed as:
```math
\begin{aligned}
0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg)
+ \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) \\
&+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2}
+ h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\
0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\
v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\
p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) .
\end{aligned}
```
Upon defining the momentum
```math
p_{n}
= h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2}
+ \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) ,
```
we can rewrite the equations of motion as
```math
\begin{aligned}
0 &= p_{n}
+ \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg)
+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\
0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\
v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\
p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\
p_{n+1}
&= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2}
+ \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) .
\end{aligned}
```
Given $(q_{n}, p_{n})$, the first four equations can be solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method, and the last equation provides an explicit update for $p_{n+1}$.
42 changes: 25 additions & 17 deletions docs/src/modules/integrators.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,16 @@
# Integrators


## Common
## Geometric Integrator

```@autodocs
Modules = [GeometricIntegrators.Integrators]
Pages = ["integrators/abstract_coefficients.jl",
"integrators/abstract_integrator.jl",
"integrators/abstract_tableau.jl",
"integrators/integrator_cache.jl",
"integrators/integrators_common.jl",
"integrators/integrators.jl"]
Pages = [
"integrators/integrator_cache.jl",
"integrators/integrator.jl",
]
```


## Initial Guesses

```@autodocs
Expand All @@ -26,7 +23,6 @@ Pages = [
]
```


## Extrapolation Methods

The extrapolation routines are exclusively used for computing
Expand All @@ -42,7 +38,6 @@ Pages = ["integrators/extrapolation/extrapolation.jl",
"integrators/extrapolation/midpoint.jl"]
```


# Euler Integrators

```@autodocs
Expand All @@ -53,8 +48,7 @@ Pages = [
]
```


## Runge-Kutta Methods
## Runge-Kutta Integrators

```@autodocs
Modules = [GeometricIntegrators.Integrators]
Expand All @@ -73,26 +67,30 @@ Pages = ["integrators/rk/abstract.jl",
]
```


## Variational Integrators

```@autodocs
Modules = [GeometricIntegrators.Integrators]
Pages = ["integrators/vi/integrators_vprk.jl"]
Pages = ["integrators/vi/vprk_methods.jl"]
Pages = [
"integrators/vi/vi_methods.jl",
"integrators/vi/position_momentum_midpoint.jl",
"integrators/vi/position_momentum_trapezoidal.jl",
"integrators/vi/vprk_integrator.jl",
"integrators/vi/vprk_methods.jl",
]
```

## Degenerate Variational Integrators

```@autodocs
Modules = [GeometricIntegrators.Integrators]
Pages = [
"integrators/dvi/methods.jl",
"integrators/dvi/integrators_dvi_a.jl",
"integrators/dvi/integrators_dvi_b.jl",
"integrators/dvi/integrators_cmdvi.jl",
"integrators/dvi/integrators_ctdvi.jl",
"integrators/dvi/integrators_dvrk.jl",
"integrators/dvi/methods.jl",
]
```

Expand All @@ -104,8 +102,18 @@ Pages = ["integrators/cgvi/integrators_cgvi.jl",
"integrators/dgvi/integrators_dgvi.jl"]
```

## Hamilton-Pontryagin Integrators

```@autodocs
Modules = [GeometricIntegrators.Integrators]
Pages = [
"integrators/hpi/hpi_methods.jl",
"integrators/hpi/hpi_midpoint.jl",
"integrators/hpi/hpi_trapezoidal.jl",
]
```

## Splitting Methods
## Splitting and Composition Methods

```@autodocs
Modules = [GeometricIntegrators.Integrators]
Expand Down
18 changes: 16 additions & 2 deletions src/Integrators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -191,9 +191,14 @@ module Integrators

export IntegratorVPRK

include("integrators/vi/vi_methods.jl")
include("integrators/vi/position_momentum_cache.jl")
include("integrators/vi/position_momentum_midpoint.jl")
include("integrators/vi/position_momentum_trapezoidal.jl")
include("integrators/vi/position_momentum_common.jl")
include("integrators/vi/vprk_methods.jl")
include("integrators/vi/integrators_vprk_cache.jl")
include("integrators/vi/integrators_vprk.jl")
include("integrators/vi/vprk_cache.jl")
include("integrators/vi/vprk_integrator.jl")


# export IntegratorCGVI, IntegratorDGVI, IntegratorDGVIEXP,
Expand All @@ -219,6 +224,15 @@ module Integrators
include("integrators/dvi/integrators_dvrk.jl")


export HPImidpointIntegrator, HPItrapezoidalIntegrator

include("integrators/hpi/hpi_methods.jl")
include("integrators/hpi/hpi_cache.jl")
include("integrators/hpi/hpi_midpoint.jl")
include("integrators/hpi/hpi_trapezoidal.jl")
include("integrators/hpi/hpi_common.jl")


export NoProjection, projection

include("projections/methods.jl")
Expand Down
6 changes: 1 addition & 5 deletions src/integrators/dvi/integrators_cmdvi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ Degenerate variational integrator cache.
* `Θ`: implicit function evaluated on solution
* `f`: vector field of implicit function
"""
struct IntegratorCacheCMDVI{DT,D} <: ODEIntegratorCache{DT,D}
struct IntegratorCacheCMDVI{DT,D} <: IODEIntegratorCache{DT,D}
x::Vector{DT}

q::Vector{DT}
Expand Down Expand Up @@ -42,10 +42,6 @@ Midpoint Degenerate Variational Integrator.
const IntegratorCMDVI{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:CMDVI}


default_solver(::CMDVI) = Newton()
default_iguess(::CMDVI) = HermiteExtrapolation()


function Base.show(io::IO, int::IntegratorCMDVI)
print(io, "\nMidpoint Degenerate Variational Integrator with:\n")
print(io, " Timestep: $(timestep(int))\n")
Expand Down
6 changes: 1 addition & 5 deletions src/integrators/dvi/integrators_ctdvi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ Degenerate variational integrator cache.
* `Θ`: implicit function evaluated on solution
* `f`: vector field of implicit function
"""
struct IntegratorCacheCTDVI{DT,D} <: ODEIntegratorCache{DT,D}
struct IntegratorCacheCTDVI{DT,D} <: IODEIntegratorCache{DT,D}
x::Vector{DT}

q::Vector{DT}
Expand Down Expand Up @@ -47,10 +47,6 @@ Trapezoidal Degenerate Variational Integrator.
const IntegratorCTDVI{DT,TT} = GeometricIntegrator{<:Union{IODEProblem{DT,TT},LODEProblem{DT,TT}}, <:CTDVI}


default_solver(::CTDVI) = Newton()
default_iguess(::CTDVI) = HermiteExtrapolation()


function Base.show(io::IO, int::IntegratorCTDVI)
print(io, "\nTrapezoidal Degenerate Variational Integrator with:\n")
print(io, " Timestep: $(timestep(int))\n")
Expand Down
Loading

0 comments on commit 1f0da91

Please sign in to comment.