Newmark beta method

Project.toml

@@ -42,6 +42,7 @@ OrdinaryDiffEqIMEXMultistep = "9f002381-b378-40b7-97a6-27a27c83f129"
 OrdinaryDiffEqLinear = "521117fe-8c41-49f8-b3b6-30780b3f0fb5"
 OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
 OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
+OrdinaryDiffEqNewmark = "d204908a-63b9-11ef-18f5-cfec7123a93b"
 OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 OrdinaryDiffEqNordsieck = "c9986a66-5c92-4813-8696-a7ec84c806c8"
 OrdinaryDiffEqPDIRK = "5dd0a6cf-3d4b-4314-aa06-06d4e299bc89"
@@ -121,8 +122,8 @@ OrdinaryDiffEqQPRK = "1"
 OrdinaryDiffEqRKN = "1"
 OrdinaryDiffEqRosenbrock = "1"
 OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqSSPRK = "1"
+OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqStabilizedRK = "1"
 OrdinaryDiffEqSymplecticRK = "1"
 OrdinaryDiffEqTsit5 = "1"

lib/OrdinaryDiffEqCore/src/algorithms.jl

@@ -46,6 +46,14 @@ abstract type OrdinaryDiffEqAdaptivePartitionedAlgorithm <: OrdinaryDiffEqAdapti
 const PartitionedAlgorithm = Union{OrdinaryDiffEqPartitionedAlgorithm,
     OrdinaryDiffEqAdaptivePartitionedAlgorithm}
 
+# Second order ODE Specific Algorithms
+abstract type OrdinaryDiffEqImplicitSecondOrderAlgorithm{CS, AD, FDT, ST, CJ} <:
+    OrdinaryDiffEqImplicitAlgorithm{CS, AD, FDT, ST, CJ} end
+abstract type OrdinaryDiffEqAdaptiveImplicitSecondOrderAlgorithm{CS, AD, FDT, ST, CJ} <:
+    OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS, AD, FDT, ST, CJ} end
+const ImplicitSecondOrderAlgorithm = Union{OrdinaryDiffEqImplicitSecondOrderAlgorithm,
+    OrdinaryDiffEqAdaptiveImplicitSecondOrderAlgorithm}
+
 function DiffEqBase.remake(thing::OrdinaryDiffEqAlgorithm; kwargs...)
     T = SciMLBase.remaker_of(thing)
     T(; SciMLBase.struct_as_namedtuple(thing)..., kwargs...)

lib/OrdinaryDiffEqNewmark/Project.toml

@@ -0,0 +1,46 @@
+name = "OrdinaryDiffEqNewmark"
+uuid = "d204908a-63b9-11ef-18f5-cfec7123a93b"
+authors = ["Dennis Ogiermann <[email protected]>"]
+version = "0.0.1"
+
+[deps]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
+
+[compat]
+DiffEqBase = "6.152.2"
+DiffEqDevTools = "2.44.4"
+FastBroadcast = "0.3.5"
+LinearAlgebra = "<0.0.1, 1"
+MacroTools = "0.5.13"
+MuladdMacro = "0.2.4"
+OrdinaryDiffEqCore = "1.1"
+OrdinaryDiffEqDifferentiation = "<0.0.1, 1"
+OrdinaryDiffEqNonlinearSolve = "<0.0.1, 1"
+Random = "<0.0.1, 1"
+RecursiveArrayTools = "3.27.0"
+Reexport = "1.2.2"
+SafeTestsets = "0.1.0"
+SciMLBase = "2.48.1"
+Test = "<0.0.1, 1"
+TruncatedStacktraces = "1.4.0"
+julia = "1.10"
+
+[extras]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"]

lib/OrdinaryDiffEqNewmark/src/OrdinaryDiffEqNewmark.jl

@@ -0,0 +1,40 @@
+module OrdinaryDiffEqNewmark
+
+import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
+                           initialize!, perform_step!, @unpack, unwrap_alg,
+                           calculate_residuals, alg_extrapolates,
+                           OrdinaryDiffEqAlgorithm,
+                           OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
+                           OrdinaryDiffEqNewtonAdaptiveAlgorithm,
+                           OrdinaryDiffEqNewtonAlgorithm,
+                           OrdinaryDiffEqImplicitSecondOrderAlgorithm,
+                           OrdinaryDiffEqAdaptiveImplicitSecondOrderAlgorithm,
+                           DEFAULT_PRECS,
+                           OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
+                           alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
+                           constvalue, _unwrap_val, _ode_interpolant,
+                           trivial_limiter!, _ode_interpolant!,
+                           isesdirk, issplit,
+                           ssp_coefficient, get_fsalfirstlast, generic_solver_docstring
+using TruncatedStacktraces, MuladdMacro, MacroTools, FastBroadcast, RecursiveArrayTools
+using SciMLBase: DynamicalODEFunction
+using LinearAlgebra: mul!, I
+import OrdinaryDiffEqCore
+
+using OrdinaryDiffEqDifferentiation: UJacobianWrapper, dolinsolve
+using OrdinaryDiffEqNonlinearSolve: du_alias_or_new, markfirststage!, build_nlsolver,
+                                    nlsolve!, nlsolvefail, isnewton, get_W, set_new_W!,
+                                    NLNewton, COEFFICIENT_MULTISTEP
+
+using Reexport
+@reexport using DiffEqBase
+
+include("algorithms.jl")
+include("alg_utils.jl")
+include("newmark_caches.jl")
+include("newmark_nlsolve.jl")
+include("newmark_perform_step.jl")
+
+export NewmarkBeta
+
+end

lib/OrdinaryDiffEqNewmark/src/alg_utils.jl

@@ -0,0 +1,5 @@
+alg_extrapolates(alg::NewmarkBeta) = true
+
+alg_order(alg::NewmarkBeta) = 1
+
+is_mass_matrix_alg(alg::NewmarkBeta) = true

lib/OrdinaryDiffEqNewmark/src/algorithms.jl

@@ -0,0 +1,53 @@
+
+"""
+    NewmarkBeta
+
+Classical Newmark-β method to solve second order ODEs, possibly in mass matrix form.
+Local truncation errors are estimated with the estimate of Zienkiewicz and Xie.
+
+## References
+
+Newmark, Nathan (1959), "A method of computation for structural dynamics",
+Journal of the Engineering Mechanics Division, 85 (EM3) (3): 67–94, doi:
+https://doi.org/10.1061/JMCEA3.0000098
+
+Zienkiewicz, O. C., and Y. M. Xie. "A simple error estimator and adaptive
+time stepping procedure for dynamic analysis." Earthquake engineering &
+structural dynamics 20.9 (1991): 871-887, doi:
+https://doi.org/10.1002/eqe.4290200907
+"""
+struct NewmarkBeta{PT, F, F2, P, CS, AD, FDT, ST, CJ} <:
+    OrdinaryDiffEqAdaptiveImplicitSecondOrderAlgorithm{CS, AD, FDT, ST, CJ}
+    β::PT
+    γ::PT
+    linsolve::F
+    nlsolve::F2
+    precs::P
+end
+
+function NewmarkBeta(β, γ; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
+    concrete_jac = nothing, diff_type = Val{:forward},
+    linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+    extrapolant = :linear)
+    NewmarkBeta{
+        typeof(β), typeof(linsolve), typeof(nlsolve), typeof(precs),
+        _unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
+        β, γ,
+        linsolve,
+        nlsolve,
+        precs)
+end
+
+# Needed for remake
+function NewmarkBeta(; β=0.25, γ=0.5, chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
+    concrete_jac = nothing, diff_type = Val{:forward},
+    linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+    extrapolant = :linear)
+    NewmarkBeta{
+        typeof(β), typeof(linsolve), typeof(nlsolve), typeof(precs),
+        _unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}(
+        β, γ,
+        linsolve,
+        nlsolve,
+        precs)
+end

lib/OrdinaryDiffEqNewmark/src/newmark_caches.jl

@@ -0,0 +1,34 @@
+@cache struct NewmarkBetaCache{uType, rateType, parameterType, N} <: OrdinaryDiffEqMutableCache
+    u::uType # Current solution
+    uprev::uType # Previous solution
+    upred::uType # Predictor solution
+    fsalfirst::rateType
+    β::parameterType # newmark parameter 1
+    γ::parameterType # newmark parameter 2
+    nlsolver::N # Inner solver
+    tmp::uType # temporary, because it is required.
+end
+
+function alg_cache(alg::NewmarkBeta, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+
+    β = alg.β
+    γ = alg.γ
+    upred = zero(u)
+    fsalfirst = zero(rate_prototype)
+
+    @assert 0.0 ≤ β ≤ 0.5
+    @assert 0.0 ≤ γ ≤ 1.0
+
+    c = 1.0
+    γ̂ = NaN # FIXME
+    nlsolver = build_nlsolver(alg, u.x[1], uprev.x[1], p, t, dt, f.f1, rate_prototype.x[1], uEltypeNoUnits,
+        uBottomEltypeNoUnits, tTypeNoUnits, γ̂, c, Val(true))
+
+    tmp = zero(u)
+    NewmarkBetaCache(u, uprev, upred, fsalfirst, β, γ, nlsolver, tmp)
+end
+
+get_fsalfirstlast(cache::NewmarkBetaCache, u) = (cache.fsalfirst, u) # FIXME use fsal

lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl

@@ -0,0 +1,106 @@
+function initialize!(integrator, cache::NewmarkBetaCache)
+    duprev, uprev = integrator.uprev.x
+    integrator.f(cache.fsalfirst, integrator.uprev, integrator.p, integrator.t)
+    integrator.stats.nf += 1
+    integrator.fsalfirst = cache.fsalfirst
+    integrator.fsallast  = cache.fsalfirst
+    # integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)
+    return
+end
+
+@muladd function perform_step!(integrator, cache::NewmarkBetaCache, repeat_step = false)
+    @unpack t, dt, f, p = integrator
+    @unpack upred, β, γ, nlsolver = cache
+
+    # This is derived from the idea stated in Nonlinear Finite Elements by Peter Wriggers, Ch 6.1.2 .
+    #
+    # Let us introduce the notation v = u' and a = u'' = v' such that we write the ODE problem as Ma = f(v,u,t).
+    # For the time discretization we assume that:
+    #   uₙ₊₁ = uₙ + Δtₙ vₙ + Δtₙ²/2 aₙ₊ₐ₁
+    #   vₙ₊₁ = vₙ + Δtₙ aₙ₊ₐ₂
+    # with a₁ = 1-2β and a₂ = 1-γ, such that
+    #   uₙ₊₁ = uₙ + Δtₙ vₙ + Δtₙ²/2 [(1-2β)aₙ + 2βaₙ₊₁]
+    #   vₙ₊₁ = vₙ + Δtₙ [(1-γ)aₙ + γaₙ₊₁]
+    # 
+    # This allows us to reduce the implicit discretization to have only aₙ₊₁ as the unknown:
+    #   Maₙ₊₁ = f(vₙ₊₁(aₙ₊₁), uₙ₊₁(aₙ₊₁), tₙ₊₁) 
+    #         = f(vₙ + Δtₙ [(1-γ)aₙ + γaₙ₊₁], uₙ + Δtₙ vₙ + Δtₙ²/2 [(1-2β)aₙ + 2βaₙ₊₁], tₙ₊₁)
+    # Such that we have to solve the nonlinear problem
+    #   Maₙ₊₁ - f(vₙ₊₁(aₙ₊₁), uₙ₊₁(aₙ₊₁), tₙ₊₁)  = 0
+    # for aₙ₊₁'' in each time step.
+
+    # For the Newton method the linearization becomes
+    #   M - (dₐuₙ₊₁ ∂fᵤ + dₐvₙ₊₁ ∂fᵥ) = 0
+    #   M - (Δtₙ²β  ∂fᵤ +  Δtₙγ  ∂fᵥ) = 0
+
+    M = f.mass_matrix
+
+    # Evaluate predictor
+    aₙ     = integrator.fsalfirst.x[1]
+    vₙ, uₙ = integrator.uprev.x
+
+    # _tmp = mass_matrix * @.. broadcast=false (α₁ * uprev+α₂ * uprev2)
+    # nlsolver.tmp = @.. broadcast=false _tmp/(dt * β₀)
+
+    # Note, we switch to notation closer to the SciML implemenation now. Needs to be double checked, also to be consistent with the formulation above
+    # nlsolve!(...) solves for
+    #   dt⋅f(innertmp + γ̂⋅z, p, t + c⋅dt) + outertmp = z
+    # So we rewrite the problem
+    #     u(tₙ₊₁)'' - f₁(ũ(tₙ₊₁) + u(tₙ₊₁)'' 2β Δtₙ², ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0
+    #   z = Δtₙ u(tₙ₊₁)'':
+    #     z         - Δtₙ f₁(ũ(tₙ₊₁) +         z 2β Δtₙ, ũ(tₙ₊₁)' +         z γ,t) = 0
+    #                 Δtₙ f₁(ũ(tₙ₊₁) +         z 2β Δtₙ, ũ(tₙ₊₁)' +         z γ,t) = z
+    #   γ̂ = [γ, 2β Δtₙ]:
+    #                 Δtₙ f₁(ũ(tₙ₊₁) +         z γ̂₂    , ũ(tₙ₊₁)' +         z γ̂₁   ,t) = z
+    #   innertmp = [ũ(tₙ₊₁)', ũ(tₙ₊₁)]:
+    #                 Δtₙ f₁(innertmp₂ +         z 2β Δtₙ², innertmp₁ +         z γ Δtₙ,t) = z
+    # Note: innertmp = nlsolve.tmp
+    # nlsolver.γ = ???
+    # nlsolver.tmp .= vₙ # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full
+    # nlsolver.z .= aₙ
+    # aₙ₊₁ = nlsolve!(nlsolver, integrator, cache, repeat_step) / dt
+    # nlsolvefail(nlsolver) && return
+
+    # Manually unrolled to see what needs to go where
+    aₙ₊₁ = copy(aₙ) # acceleration term
+    atmp = copy(aₙ)
+    J = zeros(length(aₙ), length(aₙ))
+    for i in 1:10 # = max iter - Newton loop for eq [1] above
+        uₙ₊₁ = uₙ + dt * vₙ + dt^2/2 * ((1-2β)*aₙ + 2β*aₙ₊₁)
+        vₙ₊₁ = vₙ + dt * ((1-γ)*aₙ + γ*aₙ₊₁)
+        # Compute residual
+        f.f1(atmp, vₙ₊₁, uₙ₊₁, p, t)
+        integrator.stats.nf += 1
+        residual = M*(aₙ₊₁ - atmp)
+        # Compute jacobian
+        f.jac(J, vₙ₊₁, uₙ₊₁, (γ*dt, β*dt*dt), p, t)
+        # Solve for increment
+        Δaₙ₊₁ = (M-J) \ residual
+        aₙ₊₁ .-= Δaₙ₊₁ # Looks like I messed up the signs somewhere :')
+        increment_norm = integrator.opts.internalnorm(Δaₙ₊₁, t)
+        increment_norm < 1e-4 && break
+        i == 10 && error("Newton diverged. ||Δaₙ₊₁||=$increment_norm")
+    end
+
+    u = ArrayPartition(
+        vₙ + dt * ((1-γ)*aₙ + γ*aₙ₊₁),
+        uₙ + dt * vₙ + dt^2/2 * ((1-2β)*aₙ + 2β*aₙ₊₁),
+    )
+
+    f(integrator.fsallast, u, p, t + dt)
+    integrator.stats.nf += 1
+    integrator.u = u
+
+    #
+    if integrator.opts.adaptive
+        if integrator.success_iter == 0
+            integrator.EEst = one(integrator.EEst)
+        else
+            # Zienkiewicz and Xie (1991) Eq. 21
+            δaₙ₊₁ = (integrator.fsallast.x[1] - aₙ₊₁)
+            integrator.EEst = dt*dt/2 * (2*β - 1/3) * integrator.opts.internalnorm(δaₙ₊₁, t)
+        end
+    end
+
+    return
+end

lib/OrdinaryDiffEqNewmark/test/runtests.jl

@@ -0,0 +1,71 @@
+using OrdinaryDiffEqNewmark, Test, RecursiveArrayTools, DiffEqDevTools, Statistics
+
+# Newmark methods with harmonic oscillator
+@testset "Harmonic Oscillator" begin
+    u0 = fill(0.0, 2)
+    v0 = ones(2)
+    function f1_harmonic!(dv, v, u, p, t)
+        dv .= -u
+    end
+    function harmonic_jac(J, v, u, weights, p, t)
+        J[1,1] = weights[1] * (0.0) + weights[2] * (-1.0)
+        J[1,2] = weights[1] * (0.0) + weights[2] * ( 0.0)
+        J[2,2] = weights[1] * (0.0) + weights[2] * (-1.0)
+        J[2,1] = weights[1] * (0.0) + weights[2] * ( 0.0)
+    end
+    function f2_harmonic!(du, v, u, p, t)
+        du .= v
+    end
+    function harmonic_analytic(y0, p, x)
+        v0, u0 = y0.x
+        ArrayPartition(-u0 * sin(x) + v0 * cos(x), u0 * cos(x) + v0 * sin(x))
+    end
+
+    ff_harmonic! = DynamicalODEFunction(f1_harmonic!, f2_harmonic!; jac=harmonic_jac, analytic = harmonic_analytic)
+    prob = DynamicalODEProblem(ff_harmonic!, v0, u0, (0.0, 5.0))
+    dts = 1.0 ./ 2.0 .^ (5:-1:0)
+
+    sim = test_convergence(dts, prob, NewmarkBeta(), dense_errors = true)
+    @test sim.𝒪est[:l2]≈2 rtol=1e-1
+end
+
+# Newmark methods with damped oscillator
+@testset "Damped Oscillator" begin
+    function damped_oscillator!(a, v, u, p, t)
+        @. a = -u - 0.5 * v
+        return nothing
+    end
+    function damped_jac(J, v, u, weights, p, t)
+        J[1,1] = weights[1] * (-0.5) + weights[2] * (-1.0)
+    end
+    function damped_oscillator_analytic(du0_u0, p, t)
+        ArrayPartition(
+            [
+                exp(-t / 4) / 15 * (15 * du0_u0[1] * cos(sqrt(15) * t / 4) -
+                sqrt(15) * (du0_u0[1] + 4 * du0_u0[2]) * sin(sqrt(15) * t / 4))
+            ], # du
+            [
+                exp(-t / 4) / 15 * (15 * du0_u0[2] * cos(sqrt(15) * t / 4) +
+                sqrt(15) * (4 * du0_u0[1] + du0_u0[2]) * sin(sqrt(15) * t / 4))
+            ]
+        )
+    end
+    ff_harmonic_damped! = DynamicalODEFunction(
+        damped_oscillator!,
+        (du, v, u, p, t) -> du = v;
+        jac=damped_jac,
+        analytic = damped_oscillator_analytic
+    )
+
+    prob = DynamicalODEProblem(ff_harmonic_damped!, [0.0], [1.0], (0.0, 10.0))
+    dts = 1.0 ./ 2.0 .^ (5:-1:0)
+
+    sim = test_convergence(dts, prob, NewmarkBeta(), dense_errors = true)
+    @test sim.𝒪est[:l2]≈2 rtol=1e-1
+
+    # TODO
+    # function damped_oscillator(v, u, p, t)
+    #     return -u - 0.5 * v
+    # end
+    # ...
+end