From 6dbdcd114eaa6b3b06a86ebb8c07716bbc992937 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 18 Oct 2023 18:31:59 +0200 Subject: [PATCH 01/11] first steps. --- src/OrdinaryDiffEq.jl | 2 ++ src/algorithms.jl | 20 +++++++++++++++++ src/caches/newmark_caches.jl | 23 +++++++++++++++++++ src/perform_step/newmark_perform_step.jl | 28 ++++++++++++++++++++++++ 4 files changed, 73 insertions(+) create mode 100644 src/caches/newmark_caches.jl create mode 100644 src/perform_step/newmark_perform_step.jl diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 4cefb0aec1..3ab91ef97d 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -399,6 +399,8 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, export SplitEuler +export Newmark + export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7 diff --git a/src/algorithms.jl b/src/algorithms.jl index a8959e0606..9fd62b6e1f 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -1210,6 +1210,26 @@ pages={113753} """ struct ERKN7 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end +""" + Newmark + +Classical Newmark Beta method to solve second order ODEs, possibly in mass matrix form. + +Fixed time step only. + +## References + +## 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 +""" +struct Newmark{PT} <: OrdinaryDiffEqPartitionedAlgorithm + β::PT + γ::PT +end + ################################################################################ # Adams Bashforth and Adams moulton methods diff --git a/src/caches/newmark_caches.jl b/src/caches/newmark_caches.jl new file mode 100644 index 0000000000..77c24f3dee --- /dev/null +++ b/src/caches/newmark_caches.jl @@ -0,0 +1,23 @@ +@cache struct NewmarkCache{uType, rateType, parameterType, N} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + β::parameterType + γ::parameterType + tmp::uType + nlsolver::N +end + +function alg_cache(alg::Newmark, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + + β = alg.β + γ = alg.γ + tmp = zero(u) + + nlsolver = build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, + uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true)) + + NewmarkCache(u, uprev, β, γ, tmp, nlsolver) +end diff --git a/src/perform_step/newmark_perform_step.jl b/src/perform_step/newmark_perform_step.jl new file mode 100644 index 0000000000..f504a507dd --- /dev/null +++ b/src/perform_step/newmark_perform_step.jl @@ -0,0 +1,28 @@ +function initialize!(integrator, cache::NewmarkCache) + integrator.kshortsize = 2 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + + duprev, uprev = integrator.uprev.x + ku = integrator.f(duprev, uprev, integrator.p, integrator.t) + integrator.stats.nf += 1 + integrator.stats.nf2 += 1 + integrator.fsalfirst = ku +end + +@muladd function perform_step!(integrator, cache::NewmarkCache, repeat_step = false) + @unpack t, dt, f, p = integrator + @unpack β, γ, nlsolver = cache + + mass_matrix = f.mass_matrix + + nlsolver.z = uprev + + # _tmp = mass_matrix * @.. broadcast=false (α₁ * uprev+α₂ * uprev2) + # nlsolver.tmp = @.. broadcast=false _tmp/(dt * β₀) + + # u = nlsolve!(nlsolver, integrator, cache, repeat_step) + # nlsolvefail(nlsolver) && return + + f(integrator.fsallast, u, p, t + dt) + integrator.stats.nf += 1 +end From b3634dc7794872171172c983603eed589227392d Mon Sep 17 00:00:00 2001 From: termi-official Date: Fri, 10 May 2024 20:45:48 +0200 Subject: [PATCH 02/11] Stuck state prototype --- src/OrdinaryDiffEq.jl | 5 +- src/alg_utils.jl | 3 + src/algorithms.jl | 70 +++++++++++++++++- src/caches/newmark_caches.jl | 37 ++++++---- src/misc_utils.jl | 15 ++++ src/nlsolve/newton.jl | 11 +++ src/nlsolve/nlsolve.jl | 3 + src/perform_step/newmark_perform_step.jl | 74 ++++++++++++++++--- .../partitioned_methods_tests.jl | 2 + 9 files changed, 189 insertions(+), 31 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index cc42033b57..4b7a516f56 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -164,6 +164,7 @@ include("caches/prk_caches.jl") include("caches/pdirk_caches.jl") include("caches/dae_caches.jl") include("caches/qprk_caches.jl") +include("caches/newmark_caches.jl") include("tableaus/low_order_rk_tableaus.jl") include("tableaus/high_order_rk_tableaus.jl") @@ -213,6 +214,8 @@ include("perform_step/pdirk_perform_step.jl") include("perform_step/dae_perform_step.jl") include("perform_step/qprk_perform_step.jl") +include("perform_step/newmark_perform_step.jl") + include("dense/generic_dense.jl") include("dense/interpolants.jl") include("dense/rosenbrock_interpolants.jl") @@ -404,7 +407,7 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, export SplitEuler -export Newmark +export NewmarkBeta export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, diff --git a/src/alg_utils.jl b/src/alg_utils.jl index b4acdaa7ec..7043b63e70 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -495,6 +495,7 @@ alg_order(alg::ERKN4) = 4 alg_order(alg::ERKN5) = 5 alg_order(alg::ERKN7) = 7 alg_order(alg::RKN4) = 4 +alg_order(alg::NewmarkBeta) = 1 # or at least 2, depending on the parameters. alg_order(alg::Midpoint) = 2 @@ -1064,3 +1065,5 @@ is_mass_matrix_alg(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false is_mass_matrix_alg(alg::CompositeAlgorithm) = all(is_mass_matrix_alg, alg.algs) is_mass_matrix_alg(alg::RosenbrockAlgorithm) = true is_mass_matrix_alg(alg::NewtonAlgorithm) = !isesdirk(alg) +is_mass_matrix_alg(alg::NewmarkBeta) = true + diff --git a/src/algorithms.jl b/src/algorithms.jl index 4e9d241a7a..8532665b91 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -47,8 +47,14 @@ abstract type DAEAlgorithm{CS, AD, FDT, ST, CJ} <: DiffEqBase.AbstractDAEAlgorit # Partitioned ODE Specific Algorithms abstract type OrdinaryDiffEqPartitionedAlgorithm <: OrdinaryDiffEqAlgorithm end abstract type OrdinaryDiffEqAdaptivePartitionedAlgorithm <: OrdinaryDiffEqAdaptiveAlgorithm end +abstract type OrdinaryDiffEqImplicitPartitionedAlgorithm{CS, AD, FDT, ST, CJ} <: + OrdinaryDiffEqImplicitAlgorithm{CS, AD, FDT, ST, CJ} end +abstract type OrdinaryDiffEqAdaptivePartitionedImplicitAlgorithm{CS, AD, FDT, ST, CJ} <: + OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS, AD, FDT, ST, CJ} end const PartitionedAlgorithm = Union{OrdinaryDiffEqPartitionedAlgorithm, - OrdinaryDiffEqAdaptivePartitionedAlgorithm} + OrdinaryDiffEqAdaptivePartitionedAlgorithm, + OrdinaryDiffEqImplicitPartitionedAlgorithm, + OrdinaryDiffEqAdaptivePartitionedImplicitAlgorithm} struct FunctionMap{scale_by_time} <: OrdinaryDiffEqAlgorithm end FunctionMap(; scale_by_time = false) = FunctionMap{scale_by_time}() @@ -1232,25 +1238,81 @@ year = {2024}, """ struct RKN4 <: OrdinaryDiffEqAlgorithm end - """ NewmarkBeta Classical Newmark-β method to solve second order ODEs, possibly in mass matrix form. -Fixed time step only. +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} <: OrdinaryDiffEqPartitionedAlgorithm +struct NewmarkBeta{PT, F, F2, P, CS, AD, FDT, ST, CJ} <: + OrdinaryDiffEqAdaptivePartitionedImplicitAlgorithm{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(; β=nothing, γ=nothing, 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 + +# struct CNAB2{CS, AD, F, F2, P, FDT, ST, CJ} <: +# OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} +# linsolve::F +# nlsolve::F2 +# precs::P +# extrapolant::Symbol +# end + +# function CNAB2(; 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) +# CNAB2{ +# _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), +# typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}( +# linsolve, +# nlsolve, +# precs, +# extrapolant) +# end + ################################################################################ # Adams Bashforth and Adams moulton methods diff --git a/src/caches/newmark_caches.jl b/src/caches/newmark_caches.jl index 77c24f3dee..7724c0dd51 100644 --- a/src/caches/newmark_caches.jl +++ b/src/caches/newmark_caches.jl @@ -1,23 +1,32 @@ -@cache struct NewmarkCache{uType, rateType, parameterType, N} <: OrdinaryDiffEqMutableCache - u::uType - uprev::uType - β::parameterType - γ::parameterType - tmp::uType - nlsolver::N +@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::Newmark, u, rate_prototype, ::Type{uEltypeNoUnits}, +function alg_cache(alg::NewmarkBeta, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, - ::Val{true}) {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} β = alg.β γ = alg.γ - tmp = zero(u) + upred = zero(u) + fsalfirst = zero(rate_prototype) + + @assert 0.0 ≤ β ≤ 0.5 + @assert 0.0 ≤ γ ≤ 1.0 - nlsolver = build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, - uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true)) - - NewmarkCache(u, uprev, β, γ, tmp, nlsolver) + c = 1.0 + γ̂ = ArrayPartitionNLSolveHelper(1.0,1.0) + 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 diff --git a/src/misc_utils.jl b/src/misc_utils.jl index aab109f6c6..5094b5a08b 100644 --- a/src/misc_utils.jl +++ b/src/misc_utils.jl @@ -197,3 +197,18 @@ function get_differential_vars(f, u) end end end + +""" +Helper to evaluate + f(innertmp + γ̂⋅z, p, t + c⋅dt) +in the nonlinear solver when z is an array partition and γ̂ is the scaling parameter. +Note that γ̂ is different from the γ in the Newmark-β method! +""" +struct ArrayPartitionNLSolveHelper{T} + γ₁::T + γ₂::T +end + +function Base.:*(γ::ArrayPartitionNLSolveHelper{T1}, z::ArrayPartition{T2, <: Tuple{<:Any, <:Any}}) where {T1, T2} + ArrayPartition(γ.γ₁*z.x[1], γ.γ₂*z.x[2]) +end diff --git a/src/nlsolve/newton.jl b/src/nlsolve/newton.jl index c28db2089a..87ffc7e27b 100644 --- a/src/nlsolve/newton.jl +++ b/src/nlsolve/newton.jl @@ -356,6 +356,17 @@ function compute_ustep!(ustep, tmp, γ, z, method) ustep end +# Dispatch for Newmark-type solvers +function compute_ustep!(ustep::ArrayPartition, tmp::ArrayPartition, γ::ArrayPartitionNLSolveHelper, z, method) + if method === COEFFICIENT_MULTISTEP + # ustep = z + @error "Multistep for second order ODEs does not work yet." + else + @.. ustep = ArrayPartition(tmp.x[1] + γ.γ₁ * z, tmp.x[1] + γ.γ₂ * z) + end + ustep +end + function _compute_rhs(tmp, γ, α, tstep, invγdt, method::MethodType, p, dt, f, z) mass_matrix = f.mass_matrix ustep = compute_ustep(tmp, γ, z, method) diff --git a/src/nlsolve/nlsolve.jl b/src/nlsolve/nlsolve.jl index 779ced90c9..9dabbeab6e 100644 --- a/src/nlsolve/nlsolve.jl +++ b/src/nlsolve/nlsolve.jl @@ -10,6 +10,9 @@ dt⋅f(innertmp + γ⋅z, p, t + c⋅dt) + outertmp = z ``` where `dt` is the step size and `γ` and `c` are constants, and return the solution `z`. + +Whether `innertmp` and `outertmp` is used for the evaluation is controlled by setting `nlsolver.method`. +In both cases the variable name is actually `nlsolver.tmp`. """ function nlsolve!(nlsolver::AbstractNLSolver, integrator::DiffEqBase.DEIntegrator, cache = nothing, repeat_step = false) diff --git a/src/perform_step/newmark_perform_step.jl b/src/perform_step/newmark_perform_step.jl index f504a507dd..cadd1a2625 100644 --- a/src/perform_step/newmark_perform_step.jl +++ b/src/perform_step/newmark_perform_step.jl @@ -1,28 +1,78 @@ -function initialize!(integrator, cache::NewmarkCache) - integrator.kshortsize = 2 - integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) - +function initialize!(integrator, cache::NewmarkBetaCache) duprev, uprev = integrator.uprev.x - ku = integrator.f(duprev, uprev, integrator.p, integrator.t) + integrator.f(cache.fsalfirst, integrator.uprev, integrator.p, integrator.t) + integrator.fsalfirst = cache.fsalfirst + integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst) integrator.stats.nf += 1 - integrator.stats.nf2 += 1 - integrator.fsalfirst = ku + return end -@muladd function perform_step!(integrator, cache::NewmarkCache, repeat_step = false) +@muladd function perform_step!(integrator, cache::NewmarkBetaCache, repeat_step = false) @unpack t, dt, f, p = integrator - @unpack β, γ, nlsolver = cache + @unpack upred, β, γ, nlsolver = cache + + # Given Mu'' = f(u,u',t) we need to solve the non-linear problem + # Mu(tₙ₊₁)'' - f(ũ(tₙ₊₁) + u(tₙ₊₁)'' β Δtₙ²,ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0 + # for u(tₙ₊₁)''. + # The predictors ũ(tₙ₊₁) are computed as + # ũ(tₙ₊₁) = u(tₙ) + u(tₙ)' Δtₙ + u(tₙ)'' (0.5 - β) Δtₙ² + # ũ(tₙ₊₁)' = u(tₙ)' + u(tₙ)'' (1.0 - γ) Δtₙ + # such that we can compute the solution with the correctors + # u(tₙ₊₁) = ũ(tₙ₊₁) + u(tₙ₊₁)'' β Δtₙ² + # u(tₙ₊₁)' = ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ mass_matrix = f.mass_matrix - nlsolver.z = uprev + # Evaluate predictor + dduprev = integrator.fsalfirst.x[1] + duprev, uprev = integrator.uprev.x + upred_full = ArrayPartition( + duprev + dt*(1.0 - γ)*dduprev, + uprev + dt*dt*(0.5 - β)*dduprev + dt*duprev + ) # _tmp = mass_matrix * @.. broadcast=false (α₁ * uprev+α₂ * uprev2) # nlsolver.tmp = @.. broadcast=false _tmp/(dt * β₀) - # u = nlsolve!(nlsolver, integrator, cache, repeat_step) - # nlsolvefail(nlsolver) && return + # nlsolve!(...) solves for + # dt⋅f(innertmp + γ̂⋅z, p, t + c⋅dt) + outertmp = z + # So we rewrite the problem + # u(tₙ₊₁)'' - f(ũ(tₙ₊₁) + u(tₙ₊₁)'' β Δtₙ², ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0 + # z = Δtₙ u(tₙ₊₁)'': + # z - Δtₙ f(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = 0 + # Δtₙ f(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = z + # γ̂ = [γ, β Δtₙ]: + # Δtₙ f(ũ(tₙ₊₁) + z γ̂₂ , ũ(tₙ₊₁)' + z γ̂₁ ,t) = z + # innertmp = [ũ(tₙ₊₁)', ũ(tₙ₊₁)]: + # Δtₙ f(innertmp₂ + z β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z + # Note: innertmp = nlsolve.tmp + nlsolver.γ = ArrayPartitionNLSolveHelper(γ, β * dt) # = γ̂ + nlsolver.tmp .= upred_full # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full + + # Use the linear extrapolation Δtₙ u(tₙ)'' as initial guess for the nonlinear solve + nlsolver.z .= dt*dduprev + ddu = nlsolve!(nlsolver, integrator, cache, repeat_step) / dt + nlsolvefail(nlsolver) && return + + # Apply corrector + u .= ArrayPartition( + upred_full.x[1] + ddu*β*dt*dt, + upred_full.x[2] + ddu*γ*dt + ) + + # + if integrator.opts.adaptive + if integrator.success_iter == 0 + integrator.EEst = one(integrator.EEst) + else + # Zienkiewicz and Xie (1991) Eq. 21 + δddu = (ddu - dduprev) + integrator.EEst = dt*dt/2 * (2*β - 1/3) * integrator.opts.internalnorm(δddu, t) + end + end f(integrator.fsallast, u, p, t + dt) integrator.stats.nf += 1 + + return end diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 1ef7dee128..229387c7ab 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -19,6 +19,7 @@ prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) sol = solve(prob, SymplecticEuler(), dt = 1 / 2) sol_verlet = solve(prob, VelocityVerlet(), dt = 1 / 100) sol_ruth3 = solve(prob, Ruth3(), dt = 1 / 100) +sol_newmark = solve(prob, NewmarkBeta(0.5,0.25)) interp_time = 0:0.001:5 interp = sol(0.5) @@ -31,6 +32,7 @@ prob = SecondOrderODEProblem(f1_harmonic, v0, u0, (0.0, 5.0)) sol2 = solve(prob, SymplecticEuler(), dt = 1 / 2) sol2_verlet = solve(prob, VelocityVerlet(), dt = 1 / 100) sol2_ruth3 = solve(prob, Ruth3(), dt = 1 / 100) +sol2_newmark = solve(prob, NewmarkBeta(0.5,0.25)) sol2_verlet(0.1) From d7d479e8fa8776fae779abc1fbf07732ada69492 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 15 May 2024 19:07:09 +0200 Subject: [PATCH 03/11] Almost first version. Jacobian interaction missing. --- src/algorithms.jl | 2 +- src/caches/newmark_caches.jl | 2 +- src/integrators/integrator_utils.jl | 1 + src/misc_utils.jl | 6 +++ src/nlsolve/type.jl | 10 ++-- src/nlsolve/utils.jl | 11 ++-- src/perform_step/newmark_perform_step.jl | 29 ++++++----- .../partitioned_methods_tests.jl | 51 ++++++++++++++----- 8 files changed, 74 insertions(+), 38 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 8532665b91..8295cf389d 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -3328,7 +3328,7 @@ const MultistepAlgorithms = Union{IRKN3, IRKN4, ABDF2, AB3, AB4, AB5, ABM32, ABM43, ABM54} -const SplitAlgorithms = Union{CNAB2, CNLF2, IRKC, SBDF, +const SplitAlgorithms = Union{CNAB2, CNLF2, IRKC, SBDF, NewmarkBeta, KenCarp3, KenCarp4, KenCarp47, KenCarp5, KenCarp58, CFNLIRK3} #= diff --git a/src/caches/newmark_caches.jl b/src/caches/newmark_caches.jl index 7724c0dd51..f4862950cd 100644 --- a/src/caches/newmark_caches.jl +++ b/src/caches/newmark_caches.jl @@ -25,7 +25,7 @@ function alg_cache(alg::NewmarkBeta, u, rate_prototype, ::Type{uEltypeNoUnits}, c = 1.0 γ̂ = ArrayPartitionNLSolveHelper(1.0,1.0) 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)) + uBottomEltypeNoUnits, tTypeNoUnits, γ̂, c, Val(true); tmp = zero(u)) tmp = zero(u) NewmarkBetaCache(u, uprev, upred, fsalfirst, β, γ, nlsolver, tmp) diff --git a/src/integrators/integrator_utils.jl b/src/integrators/integrator_utils.jl index 67d9ee0593..aa2c0f143e 100644 --- a/src/integrators/integrator_utils.jl +++ b/src/integrators/integrator_utils.jl @@ -507,6 +507,7 @@ nlsolve_f(f, alg::DAEAlgorithm) = f function nlsolve_f(integrator::ODEIntegrator) nlsolve_f(integrator.f, unwrap_alg(integrator, true)) end +nlsolve_f(f::DynamicalODEFunction, alg::NewmarkBeta) = f.f1 # FIXME function (integrator::ODEIntegrator)(t, ::Type{deriv} = Val{0}; idxs = nothing) where {deriv} diff --git a/src/misc_utils.jl b/src/misc_utils.jl index 5094b5a08b..c4954f363f 100644 --- a/src/misc_utils.jl +++ b/src/misc_utils.jl @@ -212,3 +212,9 @@ end function Base.:*(γ::ArrayPartitionNLSolveHelper{T1}, z::ArrayPartition{T2, <: Tuple{<:Any, <:Any}}) where {T1, T2} ArrayPartition(γ.γ₁*z.x[1], γ.γ₂*z.x[2]) end + +Base.:*(γ::ArrayPartitionNLSolveHelper{T}, scalar::T) where T = scalar*γ + +function Base.:*(scalar::T, γ::ArrayPartitionNLSolveHelper{T}) where T + ArrayPartitionNLSolveHelper(scalar*γ.γ₁, scalar*γ.γ₂) +end diff --git a/src/nlsolve/type.jl b/src/nlsolve/type.jl index bd70bef69e..5a8aae8fc7 100644 --- a/src/nlsolve/type.jl +++ b/src/nlsolve/type.jl @@ -89,11 +89,11 @@ end abstract type AbstractNLSolver{algType, iip} end -mutable struct NLSolver{algType, iip, uType, gamType, tmpType, tType, +mutable struct NLSolver{algType, iip, uType, gamType, tmpType, tmp2Type, tType, C <: AbstractNLSolverCache} <: AbstractNLSolver{algType, iip} z::uType - tmp::uType # DIRK and multistep methods only use tmp - tmp2::tmpType # for GLM if neccssary + tmp::tmpType # DIRK and multistep methods only use tmp + tmp2::tmp2Type # for GLM if neccssary ztmp::uType γ::gamType c::tType @@ -114,7 +114,7 @@ end function NLSolver{iip, tType}(z, tmp, ztmp, γ, c, α, alg, κ, fast_convergence_cutoff, ηold, iter, maxiters, status, cache, method = DIRK, tmp2 = nothing, nfails::Int = 0) where {iip, tType} - NLSolver{typeof(alg), iip, typeof(z), typeof(γ), typeof(tmp2), tType, typeof(cache)}(z, + NLSolver{typeof(alg), iip, typeof(z), typeof(γ), typeof(tmp), typeof(tmp2), tType, typeof(cache)}(z, tmp, tmp2, ztmp, @@ -157,7 +157,7 @@ mutable struct NLNewtonCache{ firststage::Bool firstcall::Bool W_γdt::tType - du1::uType + du1::rateType uf::ufType jac_config::jcType linsolve::lsType diff --git a/src/nlsolve/utils.jl b/src/nlsolve/utils.jl index 2e4cfb849c..01dce276f4 100644 --- a/src/nlsolve/utils.jl +++ b/src/nlsolve/utils.jl @@ -132,19 +132,19 @@ function build_nlsolver(alg, u, uprev, p, t, dt, f::F, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, γ, c, - iip) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + iip; tmp = zero(u)) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, - tTypeNoUnits, γ, c, 1, iip) + tTypeNoUnits, γ, c, 1, iip; tmp=tmp) end function build_nlsolver(alg, u, uprev, p, t, dt, f::F, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, γ, c, α, - iip) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + iip; tmp = zero(u)) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} build_nlsolver(alg, alg.nlsolve, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, - uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, iip) + uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, iip; tmp=tmp) end function build_nlsolver(alg, nlalg::Union{NLFunctional, NLAnderson, NLNewton, NonlinearSolveAlg}, @@ -152,7 +152,7 @@ function build_nlsolver(alg, nlalg::Union{NLFunctional, NLAnderson, NLNewton, No f::F, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, γ, c, α, - ::Val{true}) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, + ::Val{true}; tmp = zero(u)) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} #TODO #nlalg = DiffEqBase.handle_defaults(alg, nlalg) @@ -162,7 +162,6 @@ function build_nlsolver(alg, nlalg::Union{NLFunctional, NLAnderson, NLNewton, No # define fields of non-linear solver z = zero(u) - tmp = zero(u) ztmp = zero(u) # build cache of non-linear solver diff --git a/src/perform_step/newmark_perform_step.jl b/src/perform_step/newmark_perform_step.jl index cadd1a2625..2c34d47841 100644 --- a/src/perform_step/newmark_perform_step.jl +++ b/src/perform_step/newmark_perform_step.jl @@ -1,9 +1,10 @@ function initialize!(integrator, cache::NewmarkBetaCache) duprev, uprev = integrator.uprev.x integrator.f(cache.fsalfirst, integrator.uprev, integrator.p, integrator.t) - integrator.fsalfirst = cache.fsalfirst - integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst) integrator.stats.nf += 1 + integrator.fsalfirst = cache.fsalfirst + integrator.fsallast = cache.fsalfirst + # integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst) return end @@ -37,42 +38,44 @@ end # nlsolve!(...) solves for # dt⋅f(innertmp + γ̂⋅z, p, t + c⋅dt) + outertmp = z # So we rewrite the problem - # u(tₙ₊₁)'' - f(ũ(tₙ₊₁) + u(tₙ₊₁)'' β Δtₙ², ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0 + # u(tₙ₊₁)'' - f₁(ũ(tₙ₊₁) + u(tₙ₊₁)'' β Δtₙ², ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0 # z = Δtₙ u(tₙ₊₁)'': - # z - Δtₙ f(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = 0 - # Δtₙ f(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = z + # z - Δtₙ f₁(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = 0 + # Δtₙ f₁(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = z # γ̂ = [γ, β Δtₙ]: - # Δtₙ f(ũ(tₙ₊₁) + z γ̂₂ , ũ(tₙ₊₁)' + z γ̂₁ ,t) = z + # Δtₙ f₁(ũ(tₙ₊₁) + z γ̂₂ , ũ(tₙ₊₁)' + z γ̂₁ ,t) = z # innertmp = [ũ(tₙ₊₁)', ũ(tₙ₊₁)]: - # Δtₙ f(innertmp₂ + z β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z + # Δtₙ f₁(innertmp₂ + z β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z # Note: innertmp = nlsolve.tmp nlsolver.γ = ArrayPartitionNLSolveHelper(γ, β * dt) # = γ̂ nlsolver.tmp .= upred_full # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full # Use the linear extrapolation Δtₙ u(tₙ)'' as initial guess for the nonlinear solve - nlsolver.z .= dt*dduprev + @show nlsolver.z + @show dduprev + nlsolver.z = dt*dduprev ddu = nlsolve!(nlsolver, integrator, cache, repeat_step) / dt nlsolvefail(nlsolver) && return # Apply corrector - u .= ArrayPartition( + @. u = ArrayPartition( upred_full.x[1] + ddu*β*dt*dt, upred_full.x[2] + ddu*γ*dt ) + f(integrator.fsallast, u, p, t + dt) + integrator.stats.nf += 1 + # if integrator.opts.adaptive if integrator.success_iter == 0 integrator.EEst = one(integrator.EEst) else # Zienkiewicz and Xie (1991) Eq. 21 - δddu = (ddu - dduprev) + @. δddu = (integrator.fsallast - ddu) integrator.EEst = dt*dt/2 * (2*β - 1/3) * integrator.opts.internalnorm(δddu, t) end end - f(integrator.fsallast, u, p, t + dt) - integrator.stats.nf += 1 - return end diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 229387c7ab..b5529274b0 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -19,6 +19,11 @@ prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) sol = solve(prob, SymplecticEuler(), dt = 1 / 2) sol_verlet = solve(prob, VelocityVerlet(), dt = 1 / 100) sol_ruth3 = solve(prob, Ruth3(), dt = 1 / 100) + +harmonic_jac_f1!(J, v, u, p, t) = J[1,1] = -1.0 +ff_f1 = ODEFunction{false}(ODEFunction{false}(f1_harmonic); jac=harmonic_jac_f1!) +ff_harmonic = DynamicalODEFunction(ff_f1, f2_harmonic; analytic = harmonic_analytic) +prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) sol_newmark = solve(prob, NewmarkBeta(0.5,0.25)) interp_time = 0:0.001:5 @@ -387,26 +392,48 @@ sol = solve(prob, ERKN5(), reltol = 1e-8) sol = solve(prob, ERKN7(), reltol = 1e-8) @test length(sol.u) < 38 +# Newmark methods with damped oscillator + +ff_harmonic_damped = DynamicalODEFunction((ddu, v, u, p, t) -> ddu = -u - 0.5 * v, + (du, v, u, p, t) -> du = v, + analytic = (du0_u0, p, t) -> OrdinaryDiffEq.SciMLBase.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)) + ]) +) + +prob = DynamicalODEProblem(ff_harmonic_damped, [0.0], [1.0], (0.0, 10.0)) +sol_newmark = solve(prob, NewmarkBeta(0.5,0.25)) + +ff_harmonic_damped = DynamicalODEFunction{false}((du, u, p, t) -> -u - 0.5 * du, + (du, u, p, t) -> du, + analytic = (du0_u0, p, t) -> OrdinaryDiffEq.SciMLBase.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)) + ]) +) + # Testing generalized Runge-Kutte-Nyström methods on velocity dependend ODEs with the damped oscillator println("Out of Place") # Damped oscillator prob = ODEProblem( - DynamicalODEFunction{false}((du, u, p, t) -> -u - 0.5 * du, - (du, u, p, t) -> du, - analytic = (du0_u0, p, t) -> OrdinaryDiffEq.SciMLBase.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)) - ])), + ff_harmonic_damped, OrdinaryDiffEq.SciMLBase.ArrayPartition([0.0], [1.0]), # du0, u0 (0.0, 10.0), # tspan DiffEqBase.NullParameters(), # p - SecondOrderODEProblem{false}()) + SecondOrderODEProblem{false}() +) dts = 1.0 ./ 2.0 .^ (5:-1:0) sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true) From 1380a6c5db5008b04a56313778ad81b6ac951fd1 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 15 May 2024 19:28:25 +0200 Subject: [PATCH 04/11] Missing multiply --- src/misc_utils.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/misc_utils.jl b/src/misc_utils.jl index c4954f363f..f65accd584 100644 --- a/src/misc_utils.jl +++ b/src/misc_utils.jl @@ -213,6 +213,10 @@ function Base.:*(γ::ArrayPartitionNLSolveHelper{T1}, z::ArrayPartition{T2, <: T ArrayPartition(γ.γ₁*z.x[1], γ.γ₂*z.x[2]) end +function Base.:*(γ::ArrayPartitionNLSolveHelper{T1}, z::Vector{T2}) where {T1, T2} + ArrayPartition(γ.γ₁*z, γ.γ₂*z) +end + Base.:*(γ::ArrayPartitionNLSolveHelper{T}, scalar::T) where T = scalar*γ function Base.:*(scalar::T, γ::ArrayPartitionNLSolveHelper{T}) where T From 7b4181343155ddb1d1a8c492ba4bcb6b37bf8848 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 15 May 2024 23:03:59 +0200 Subject: [PATCH 05/11] oopsie. --- src/algorithms.jl | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 8295cf389d..560dc1943d 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -1292,27 +1292,6 @@ function NewmarkBeta(; β=nothing, γ=nothing, chunk_size = Val{0}(), autodiff = precs) end -# struct CNAB2{CS, AD, F, F2, P, FDT, ST, CJ} <: -# OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} -# linsolve::F -# nlsolve::F2 -# precs::P -# extrapolant::Symbol -# end - -# function CNAB2(; 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) -# CNAB2{ -# _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), -# typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}( -# linsolve, -# nlsolve, -# precs, -# extrapolant) -# end - ################################################################################ # Adams Bashforth and Adams moulton methods From 2e65db75a79986ea2b4b2b58c3108693b5d5bbbf Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 15 May 2024 23:57:48 +0200 Subject: [PATCH 06/11] Crashes fixed. But the step does nothing, too. --- src/derivative_utils.jl | 13 ++++++++++++ src/misc_utils.jl | 2 ++ src/nlsolve/newton.jl | 17 ++++++++++++++- src/nlsolve/utils.jl | 1 + src/perform_step/newmark_perform_step.jl | 6 ++---- .../partitioned_methods_tests.jl | 21 ++++++++++++------- 6 files changed, 47 insertions(+), 13 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 272482dfc2..02c0c290c3 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -629,6 +629,19 @@ function jacobian2W(mass_matrix::MT, dtgamma::Number, J::AbstractMatrix, return W end +# The addition here does just work for the given example. We need to split up the jacobian into its blocks. +function jacobian2W(mass_matrix, dtgamma::ArrayPartitionNLSolveHelper, + J::AbstractMatrix, W_transform::Bool) + @unpack γ₁, γ₂ = dtgamma + jacobian2W(mass_matrix,γ₁+γ₂,J,W_transform) +end + +function jacobian2W!(W::AbstractMatrix, mass_matrix, dtgamma::ArrayPartitionNLSolveHelper, + J::AbstractMatrix, W_transform::Bool) + @unpack γ₁, γ₂ = dtgamma + jacobian2W!(W,mass_matrix,γ₁+γ₂,J,W_transform) +end + function calc_W!(W, integrator, nlsolver::Union{Nothing, AbstractNLSolver}, cache, dtgamma, repeat_step, W_transform = false, newJW = nothing) @unpack t, dt, uprev, u, f, p = integrator diff --git a/src/misc_utils.jl b/src/misc_utils.jl index f65accd584..90ab8f3f3c 100644 --- a/src/misc_utils.jl +++ b/src/misc_utils.jl @@ -222,3 +222,5 @@ Base.:*(γ::ArrayPartitionNLSolveHelper{T}, scalar::T) where T = scalar*γ function Base.:*(scalar::T, γ::ArrayPartitionNLSolveHelper{T}) where T ArrayPartitionNLSolveHelper(scalar*γ.γ₁, scalar*γ.γ₂) end + +Base.inv(γ::ArrayPartitionNLSolveHelper) = 1/(γ.γ₁+γ.γ₂) diff --git a/src/nlsolve/newton.jl b/src/nlsolve/newton.jl index 87ffc7e27b..41bbcbfcc4 100644 --- a/src/nlsolve/newton.jl +++ b/src/nlsolve/newton.jl @@ -303,7 +303,7 @@ end # page 54. γdt = isdae ? α * invγdt : γ * dt - !(W_γdt ≈ γdt) && (rmul!(dz, 2 / (1 + γdt / W_γdt))) + # !(W_γdt ≈ γdt) && (rmul!(dz, 2 / (1 + γdt / W_γdt))) relax!(dz, nlsolver, integrator, f) calculate_residuals!(atmp, dz, uprev, ustep, opts.abstol, opts.reltol, @@ -423,6 +423,21 @@ function _compute_rhs!(tmp, ztmp, ustep, γ, α, tstep, k, return _vec(ztmp), ustep end +function _compute_rhs!(tmp, ztmp, ustep, γ::ArrayPartitionNLSolveHelper, α, tstep, k, + invγdt, method::MethodType, p, dt, f, z) + mass_matrix = f.mass_matrix + ustep = tmp + γ * z + f(k, ustep, p, tstep) + if mass_matrix === I + @.. ztmp=(dt * k - z) * invγdt + else + update_coefficients!(mass_matrix, ustep, p, tstep) + mul!(_vec(ztmp), mass_matrix, _vec(z)) + @.. ztmp=(dt * k - ztmp) * invγdt + end + return _vec(ztmp), ustep +end + function _compute_rhs!(tmp::Array, ztmp::Array, ustep::Array, α, tstep, k, invγdt, p, uprev, f::TF, z) where {TF<:DAEFunction} @inbounds @simd ivdep for i in eachindex(z) diff --git a/src/nlsolve/utils.jl b/src/nlsolve/utils.jl index 01dce276f4..39e77f76be 100644 --- a/src/nlsolve/utils.jl +++ b/src/nlsolve/utils.jl @@ -42,6 +42,7 @@ get_new_W!(::AbstractNLSolverCache)::Bool = true get_W(nlsolver::AbstractNLSolver) = get_W(nlsolver.cache) get_W(nlcache::Union{NLNewtonCache, NLNewtonConstantCache}) = nlcache.W +set_W_γdt!(nlsolver::AbstractNLSolver, W_γdt::ArrayPartitionNLSolveHelper) = set_W_γdt!(nlsolver, W_γdt.γ₁ + W_γdt.γ₂) set_W_γdt!(nlsolver::AbstractNLSolver, W_γdt) = set_W_γdt!(nlsolver.cache, W_γdt) function set_W_γdt!(nlcache::Union{NLNewtonCache, NLNewtonConstantCache}, W_γdt) nlcache.W_γdt = W_γdt diff --git a/src/perform_step/newmark_perform_step.jl b/src/perform_step/newmark_perform_step.jl index 2c34d47841..236a3a2207 100644 --- a/src/perform_step/newmark_perform_step.jl +++ b/src/perform_step/newmark_perform_step.jl @@ -51,14 +51,12 @@ end nlsolver.tmp .= upred_full # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full # Use the linear extrapolation Δtₙ u(tₙ)'' as initial guess for the nonlinear solve - @show nlsolver.z - @show dduprev nlsolver.z = dt*dduprev ddu = nlsolve!(nlsolver, integrator, cache, repeat_step) / dt nlsolvefail(nlsolver) && return # Apply corrector - @. u = ArrayPartition( + u = ArrayPartition( upred_full.x[1] + ddu*β*dt*dt, upred_full.x[2] + ddu*γ*dt ) @@ -72,7 +70,7 @@ end integrator.EEst = one(integrator.EEst) else # Zienkiewicz and Xie (1991) Eq. 21 - @. δddu = (integrator.fsallast - ddu) + δddu = (integrator.fsallast.x[1] - ddu) integrator.EEst = dt*dt/2 * (2*β - 1/3) * integrator.opts.internalnorm(δddu, t) end end diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index b5529274b0..cb3e53508b 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -2,15 +2,15 @@ using OrdinaryDiffEq, Test, RecursiveArrayTools, DiffEqDevTools, Statistics u0 = fill(0.0, 2) v0 = ones(2) -function f1_harmonic(dv, v, u, p, t) +function f1_harmonic(dv, u, v, p, t) dv .= -u end -function f2_harmonic(du, v, u, p, t) +function f2_harmonic(du, u, v, 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)) + u0, v0 = y0.x + ArrayPartition(u0 * cos(x) + v0 * sin(x), -u0 * sin(x) + v0 * cos(x)) end ff_harmonic = DynamicalODEFunction(f1_harmonic, f2_harmonic; analytic = harmonic_analytic) prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) @@ -20,11 +20,16 @@ sol = solve(prob, SymplecticEuler(), dt = 1 / 2) sol_verlet = solve(prob, VelocityVerlet(), dt = 1 / 100) sol_ruth3 = solve(prob, Ruth3(), dt = 1 / 100) -harmonic_jac_f1!(J, v, u, p, t) = J[1,1] = -1.0 -ff_f1 = ODEFunction{false}(ODEFunction{false}(f1_harmonic); jac=harmonic_jac_f1!) +# We need these to trick ODEFunction into accepting the split function +f1_harmonic(dv, vu::ArrayPartition, p, t) = f1_harmonic(dv, vu.x[1], vu.x[2], p, t) +f2_harmonic(du, vu::ArrayPartition, p, t) = f2_harmonic(du, vu.x[1], vu.x[2], p, t) +# We use these to trick the Newton into computing the right Jacobian +harmonic_jac_f1!(J, vu::ArrayPartition, p, t) = harmonic_jac_f1!(J, vu.x[1], vu.x[2], p, t) +harmonic_jac_f1!(J, u, v, p, t) = J[1,1] = -1.0 +ff_f1 = ODEFunction(f1_harmonic; jac=harmonic_jac_f1!) ff_harmonic = DynamicalODEFunction(ff_f1, f2_harmonic; analytic = harmonic_analytic) -prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) -sol_newmark = solve(prob, NewmarkBeta(0.5,0.25)) +prob = DynamicalODEProblem(ff_harmonic, u0, v0, (0.0, 5.0)) +sol_newmark = solve(prob, NewmarkBeta(0.5,0.25), dt = 1 / 2) interp_time = 0:0.001:5 interp = sol(0.5) From 9bb240146517edcb20ac6f89ffa0053531d29972 Mon Sep 17 00:00:00 2001 From: termi-official Date: Thu, 16 May 2024 00:28:07 +0200 Subject: [PATCH 07/11] Solves pendulum. --- src/perform_step/newmark_perform_step.jl | 6 ++++-- test/algconvergence/partitioned_methods_tests.jl | 12 ++++++------ 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/perform_step/newmark_perform_step.jl b/src/perform_step/newmark_perform_step.jl index 236a3a2207..14c0fab436 100644 --- a/src/perform_step/newmark_perform_step.jl +++ b/src/perform_step/newmark_perform_step.jl @@ -48,6 +48,7 @@ end # Δtₙ f₁(innertmp₂ + z β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z # Note: innertmp = nlsolve.tmp nlsolver.γ = ArrayPartitionNLSolveHelper(γ, β * dt) # = γ̂ + # nlsolver.γ = ArrayPartitionNLSolveHelper(0.0, β * dt) nlsolver.tmp .= upred_full # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full # Use the linear extrapolation Δtₙ u(tₙ)'' as initial guess for the nonlinear solve @@ -57,12 +58,13 @@ end # Apply corrector u = ArrayPartition( - upred_full.x[1] + ddu*β*dt*dt, - upred_full.x[2] + ddu*γ*dt + upred_full.x[1] + ddu*γ*dt, + upred_full.x[2] + ddu*β*dt*dt, ) f(integrator.fsallast, u, p, t + dt) integrator.stats.nf += 1 + integrator.u = u # if integrator.opts.adaptive diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index cb3e53508b..207cd445cd 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -2,15 +2,15 @@ using OrdinaryDiffEq, Test, RecursiveArrayTools, DiffEqDevTools, Statistics u0 = fill(0.0, 2) v0 = ones(2) -function f1_harmonic(dv, u, v, p, t) +function f1_harmonic(dv, v, u, p, t) dv .= -u end -function f2_harmonic(du, u, v, p, t) +function f2_harmonic(du, v, u, p, t) du .= v end function harmonic_analytic(y0, p, x) - u0, v0 = y0.x - ArrayPartition(u0 * cos(x) + v0 * sin(x), -u0 * sin(x) + v0 * cos(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; analytic = harmonic_analytic) prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) @@ -28,8 +28,8 @@ harmonic_jac_f1!(J, vu::ArrayPartition, p, t) = harmonic_jac_f1!(J, vu.x[1], vu. harmonic_jac_f1!(J, u, v, p, t) = J[1,1] = -1.0 ff_f1 = ODEFunction(f1_harmonic; jac=harmonic_jac_f1!) ff_harmonic = DynamicalODEFunction(ff_f1, f2_harmonic; analytic = harmonic_analytic) -prob = DynamicalODEProblem(ff_harmonic, u0, v0, (0.0, 5.0)) -sol_newmark = solve(prob, NewmarkBeta(0.5,0.25), dt = 1 / 2) +prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) +sol_newmark = solve(prob, NewmarkBeta(0.5,0.25), dt = 1 / 100, adaptive=false) interp_time = 0:0.001:5 interp = sol(0.5) From c2c7db3730a2ffdfe84391ed8ccbc47385bf12b7 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 28 Aug 2024 16:26:24 +0200 Subject: [PATCH 08/11] Hand unroll the algorithm --- lib/OrdinaryDiffEqCore/src/algorithms.jl | 8 ++ .../src/OrdinaryDiffEqNewmark.jl | 2 + lib/OrdinaryDiffEqNewmark/src/algorithms.jl | 2 +- .../src/newmark_perform_step.jl | 80 ++++++++++++------- 4 files changed, 62 insertions(+), 30 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/algorithms.jl b/lib/OrdinaryDiffEqCore/src/algorithms.jl index 01c24467a6..5d0b5c61e4 100644 --- a/lib/OrdinaryDiffEqCore/src/algorithms.jl +++ b/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...) diff --git a/lib/OrdinaryDiffEqNewmark/src/OrdinaryDiffEqNewmark.jl b/lib/OrdinaryDiffEqNewmark/src/OrdinaryDiffEqNewmark.jl index 5a48a55200..c773b4c86c 100644 --- a/lib/OrdinaryDiffEqNewmark/src/OrdinaryDiffEqNewmark.jl +++ b/lib/OrdinaryDiffEqNewmark/src/OrdinaryDiffEqNewmark.jl @@ -7,6 +7,8 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, OrdinaryDiffEqNewtonAdaptiveAlgorithm, OrdinaryDiffEqNewtonAlgorithm, + OrdinaryDiffEqImplicitSecondOrderAlgorithm, + OrdinaryDiffEqAdaptiveImplicitSecondOrderAlgorithm, DEFAULT_PRECS, OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, alg_cache, _vec, _reshape, @cache, isfsal, full_cache, diff --git a/lib/OrdinaryDiffEqNewmark/src/algorithms.jl b/lib/OrdinaryDiffEqNewmark/src/algorithms.jl index e502a9b5b4..863c03798b 100644 --- a/lib/OrdinaryDiffEqNewmark/src/algorithms.jl +++ b/lib/OrdinaryDiffEqNewmark/src/algorithms.jl @@ -17,7 +17,7 @@ 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} <: - OrdinaryDiffEqAdaptivePartitionedImplicitAlgorithm{CS, AD, FDT, ST, CJ} + OrdinaryDiffEqAdaptiveImplicitSecondOrderAlgorithm{CS, AD, FDT, ST, CJ} β::PT γ::PT linsolve::F diff --git a/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl b/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl index 14c0fab436..31261bf05a 100644 --- a/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl +++ b/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl @@ -12,50 +12,72 @@ end @unpack t, dt, f, p = integrator @unpack upred, β, γ, nlsolver = cache - # Given Mu'' = f(u,u',t) we need to solve the non-linear problem - # Mu(tₙ₊₁)'' - f(ũ(tₙ₊₁) + u(tₙ₊₁)'' β Δtₙ²,ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0 - # for u(tₙ₊₁)''. - # The predictors ũ(tₙ₊₁) are computed as - # ũ(tₙ₊₁) = u(tₙ) + u(tₙ)' Δtₙ + u(tₙ)'' (0.5 - β) Δtₙ² - # ũ(tₙ₊₁)' = u(tₙ)' + u(tₙ)'' (1.0 - γ) Δtₙ - # such that we can compute the solution with the correctors - # u(tₙ₊₁) = ũ(tₙ₊₁) + u(tₙ₊₁)'' β Δtₙ² - # u(tₙ₊₁)' = ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ + # 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(u,v,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(uₙ₊₁(aₙ₊₁), vₙ₊₁(aₙ₊₁), tₙ₊₁) + # = f(uₙ + Δtₙ vₙ + Δtₙ²/2 [(1-2β)aₙ + 2βaₙ₊₁], vₙ + Δtₙ [(1-γ)aₙ + γaₙ₊₁], tₙ₊₁) + # Such that we have to solve the nonlinear problem + # Maₙ₊₁ - f(uₙ₊₁(aₙ₊₁), vₙ₊₁(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 - mass_matrix = f.mass_matrix + M = f.mass_matrix # Evaluate predictor - dduprev = integrator.fsalfirst.x[1] - duprev, uprev = integrator.uprev.x - upred_full = ArrayPartition( - duprev + dt*(1.0 - γ)*dduprev, - uprev + dt*dt*(0.5 - β)*dduprev + dt*duprev - ) + aₙ = integrator.fsalfirst.x[1] + uₙ, vₙ = 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ₙ₊₁)'' β Δtₙ², ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0 + # u(tₙ₊₁)'' - f₁(ũ(tₙ₊₁) + u(tₙ₊₁)'' 2β Δtₙ², ũ(tₙ₊₁)' + u(tₙ₊₁)'' γ Δtₙ,t) = 0 # z = Δtₙ u(tₙ₊₁)'': - # z - Δtₙ f₁(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = 0 - # Δtₙ f₁(ũ(tₙ₊₁) + z β Δtₙ, ũ(tₙ₊₁)' + z γ,t) = z - # γ̂ = [γ, β Δ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 β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z + # Δtₙ f₁(innertmp₂ + z 2β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z # Note: innertmp = nlsolve.tmp - nlsolver.γ = ArrayPartitionNLSolveHelper(γ, β * dt) # = γ̂ - # nlsolver.γ = ArrayPartitionNLSolveHelper(0.0, β * dt) - nlsolver.tmp .= upred_full # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full - - # Use the linear extrapolation Δtₙ u(tₙ)'' as initial guess for the nonlinear solve - nlsolver.z = dt*dduprev - ddu = nlsolve!(nlsolver, integrator, cache, repeat_step) / dt - nlsolvefail(nlsolver) && return + # nlsolver.γ = ??? + # nlsolver.tmp .= vₙ # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full + # nlsolver.z .= aₙ + # ddu = 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ₙ) + for i in 1:10 # = max iter - Newton loop for eq [1] above + uₙ₊₁ = uₙ + Δtₙ * vₙ + Δtₙ^2/2 * ((1-2β)*aₙ + 2β*aₙ₊₁) + vₙ₊₁ = vₙ + Δtₙ * ((1-γ)aₙ + γaₙ₊₁) + # Compute residual + f.f1(atmp, uₙ₊₁, vₙ₊₁, p, t) + residual = M*(aₙ₊₁ - atmp) + # Compute jacobian + f.jac(J, upred_full, (β*dt*dt, γ*dt), t) + # Solve for increment + Δddu = (M-J) \ residual + ddu .+= Δddu + norm(Δddu) < 1e-4 && break + i == 10 && error("Newton diverged. Δddu=$(norm(Δddu))") + end # Apply corrector u = ArrayPartition( upred_full.x[1] + ddu*γ*dt, From 68498e93a4d869f6c140de93968b1bdde6920470 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 28 Aug 2024 16:37:42 +0200 Subject: [PATCH 09/11] Looks like someone does not know how to merge. --- lib/OrdinaryDiffEqNonlinearSolve/src/type.jl | 11 ++++++----- lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl | 11 ++++++----- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl index e4d0180af9..7a0095c3f3 100644 --- a/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl +++ b/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl @@ -68,11 +68,11 @@ end # solver -mutable struct NLSolver{algType, iip, uType, gamType, tmpType, tmp2Type, tType, - C <: AbstractNLSolverCache} <: AbstractNLSolver{algType, iip} +mutable struct NLSolver{algType, iip, uType, gamType, tmpType, tType, + C <: AbstractNLSolverCache, E} <: AbstractNLSolver{algType, iip} z::uType - tmp::tmpType # DIRK and multistep methods only use tmp - tmp2::tmp2Type # for GLM if neccssary + tmp::uType # DIRK and multistep methods only use tmp + tmp2::tmpType # for GLM if neccssary ztmp::uType γ::gamType c::tType @@ -94,7 +94,8 @@ end function NLSolver{iip, tType}(z, tmp, ztmp, γ, c, α, alg, κ, fast_convergence_cutoff, ηold, iter, maxiters, status, cache, method = DIRK, tmp2 = nothing, nfails::Int = 0) where {iip, tType} - NLSolver{typeof(alg), iip, typeof(z), typeof(γ), typeof(tmp), typeof(tmp2), tType, typeof(cache)}(z, + NLSolver{typeof(alg), iip, typeof(z), typeof(γ), typeof(tmp), tType, typeof(cache), RT}( + z, tmp, tmp2, ztmp, diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl index a44dcd48d1..3bed6a453a 100644 --- a/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl +++ b/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl @@ -130,19 +130,19 @@ function build_nlsolver(alg, u, uprev, p, t, dt, f::F, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, γ, c, - iip; tmp = zero(u)) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + iip) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, - tTypeNoUnits, γ, c, 1, iip; tmp) + tTypeNoUnits, γ, c, 1, iip) end function build_nlsolver(alg, u, uprev, p, t, dt, f::F, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, γ, c, α, - iip; tmp = zero(u)) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + iip) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} build_nlsolver(alg, alg.nlsolve, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, - uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, iip; tmp) + uBottomEltypeNoUnits, tTypeNoUnits, γ, c, α, iip) end function build_nlsolver( @@ -151,7 +151,7 @@ function build_nlsolver( f::F, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, γ, c, α, - ::Val{true}; tmp = zero(u)) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, + ::Val{true}) where {F, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} #TODO #nlalg = DiffEqBase.handle_defaults(alg, nlalg) @@ -161,6 +161,7 @@ function build_nlsolver( # define fields of non-linear solver z = zero(u) + tmp = zero(u) ztmp = zero(u) # build cache of non-linear solver From 302903da5830a256d38864ac721baaa7525a3037 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 28 Aug 2024 16:39:35 +0200 Subject: [PATCH 10/11] Again. --- lib/OrdinaryDiffEqNonlinearSolve/src/type.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl index 7a0095c3f3..b6a334ca1f 100644 --- a/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl +++ b/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl @@ -94,7 +94,8 @@ end function NLSolver{iip, tType}(z, tmp, ztmp, γ, c, α, alg, κ, fast_convergence_cutoff, ηold, iter, maxiters, status, cache, method = DIRK, tmp2 = nothing, nfails::Int = 0) where {iip, tType} - NLSolver{typeof(alg), iip, typeof(z), typeof(γ), typeof(tmp), tType, typeof(cache), RT}( + RT = real(eltype(z)) + NLSolver{typeof(alg), iip, typeof(z), typeof(γ), typeof(tmp2), tType, typeof(cache), RT}( z, tmp, tmp2, From b07101ba7c9ba10fb4547c11408c22d774584df3 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 28 Aug 2024 17:19:36 +0200 Subject: [PATCH 11/11] First functional version --- lib/OrdinaryDiffEqNewmark/src/algorithms.jl | 2 +- .../src/newmark_caches.jl | 6 +- .../src/newmark_perform_step.jl | 41 +++++------ lib/OrdinaryDiffEqNewmark/test/runtests.jl | 68 ++++++++----------- 4 files changed, 56 insertions(+), 61 deletions(-) diff --git a/lib/OrdinaryDiffEqNewmark/src/algorithms.jl b/lib/OrdinaryDiffEqNewmark/src/algorithms.jl index 863c03798b..9a22a70758 100644 --- a/lib/OrdinaryDiffEqNewmark/src/algorithms.jl +++ b/lib/OrdinaryDiffEqNewmark/src/algorithms.jl @@ -39,7 +39,7 @@ function NewmarkBeta(β, γ; chunk_size = Val{0}(), autodiff = Val{true}(), stan end # Needed for remake -function NewmarkBeta(; β=nothing, γ=nothing, chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), +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) diff --git a/lib/OrdinaryDiffEqNewmark/src/newmark_caches.jl b/lib/OrdinaryDiffEqNewmark/src/newmark_caches.jl index 1a12103a67..d536ce7dda 100644 --- a/lib/OrdinaryDiffEqNewmark/src/newmark_caches.jl +++ b/lib/OrdinaryDiffEqNewmark/src/newmark_caches.jl @@ -23,10 +23,12 @@ function alg_cache(alg::NewmarkBeta, u, rate_prototype, ::Type{uEltypeNoUnits}, @assert 0.0 ≤ γ ≤ 1.0 c = 1.0 - γ̂ = ArrayPartitionNLSolveHelper(1.0,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)) + 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 diff --git a/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl b/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl index 31261bf05a..d01f736803 100644 --- a/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl +++ b/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl @@ -14,7 +14,7 @@ end # 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(u,v,t). + # 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ₙ₊ₐ₂ @@ -23,10 +23,10 @@ end # vₙ₊₁ = vₙ + Δtₙ [(1-γ)aₙ + γaₙ₊₁] # # This allows us to reduce the implicit discretization to have only aₙ₊₁ as the unknown: - # Maₙ₊₁ = f(uₙ₊₁(aₙ₊₁), vₙ₊₁(aₙ₊₁), tₙ₊₁) - # = f(uₙ + Δtₙ vₙ + Δtₙ²/2 [(1-2β)aₙ + 2βaₙ₊₁], vₙ + Δtₙ [(1-γ)aₙ + γaₙ₊₁], tₙ₊₁) + # 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(uₙ₊₁(aₙ₊₁), vₙ₊₁(aₙ₊₁), tₙ₊₁) = 0 + # Maₙ₊₁ - f(vₙ₊₁(aₙ₊₁), uₙ₊₁(aₙ₊₁), tₙ₊₁) = 0 # for aₙ₊₁'' in each time step. # For the Newton method the linearization becomes @@ -37,7 +37,7 @@ end # Evaluate predictor aₙ = integrator.fsalfirst.x[1] - uₙ, vₙ = integrator.uprev.x + vₙ, uₙ = integrator.uprev.x # _tmp = mass_matrix * @.. broadcast=false (α₁ * uprev+α₂ * uprev2) # nlsolver.tmp = @.. broadcast=false _tmp/(dt * β₀) @@ -58,30 +58,33 @@ end # nlsolver.γ = ??? # nlsolver.tmp .= vₙ # TODO check f tmp is potentially modified and if not elimiate the allocation of upred_full # nlsolver.z .= aₙ - # ddu = nlsolve!(nlsolver, integrator, cache, repeat_step) / dt + # 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ₙ + Δtₙ * vₙ + Δtₙ^2/2 * ((1-2β)*aₙ + 2β*aₙ₊₁) - vₙ₊₁ = vₙ + Δtₙ * ((1-γ)aₙ + γaₙ₊₁) + uₙ₊₁ = uₙ + dt * vₙ + dt^2/2 * ((1-2β)*aₙ + 2β*aₙ₊₁) + vₙ₊₁ = vₙ + dt * ((1-γ)*aₙ + γ*aₙ₊₁) # Compute residual - f.f1(atmp, uₙ₊₁, vₙ₊₁, p, t) + f.f1(atmp, vₙ₊₁, uₙ₊₁, p, t) + integrator.stats.nf += 1 residual = M*(aₙ₊₁ - atmp) # Compute jacobian - f.jac(J, upred_full, (β*dt*dt, γ*dt), t) + f.jac(J, vₙ₊₁, uₙ₊₁, (γ*dt, β*dt*dt), p, t) # Solve for increment - Δddu = (M-J) \ residual - ddu .+= Δddu - norm(Δddu) < 1e-4 && break - i == 10 && error("Newton diverged. Δddu=$(norm(Δddu))") + Δ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 - # Apply corrector + u = ArrayPartition( - upred_full.x[1] + ddu*γ*dt, - upred_full.x[2] + ddu*β*dt*dt, + vₙ + dt * ((1-γ)*aₙ + γ*aₙ₊₁), + uₙ + dt * vₙ + dt^2/2 * ((1-2β)*aₙ + 2β*aₙ₊₁), ) f(integrator.fsallast, u, p, t + dt) @@ -94,8 +97,8 @@ end integrator.EEst = one(integrator.EEst) else # Zienkiewicz and Xie (1991) Eq. 21 - δddu = (integrator.fsallast.x[1] - ddu) - integrator.EEst = dt*dt/2 * (2*β - 1/3) * integrator.opts.internalnorm(δddu, t) + δaₙ₊₁ = (integrator.fsallast.x[1] - aₙ₊₁) + integrator.EEst = dt*dt/2 * (2*β - 1/3) * integrator.opts.internalnorm(δaₙ₊₁, t) end end diff --git a/lib/OrdinaryDiffEqNewmark/test/runtests.jl b/lib/OrdinaryDiffEqNewmark/test/runtests.jl index 2c0bd54468..4b90eb2815 100644 --- a/lib/OrdinaryDiffEqNewmark/test/runtests.jl +++ b/lib/OrdinaryDiffEqNewmark/test/runtests.jl @@ -7,6 +7,12 @@ using OrdinaryDiffEqNewmark, Test, RecursiveArrayTools, DiffEqDevTools, Statisti 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 @@ -15,25 +21,25 @@ using OrdinaryDiffEqNewmark, Test, RecursiveArrayTools, DiffEqDevTools, Statisti ArrayPartition(-u0 * sin(x) + v0 * cos(x), u0 * cos(x) + v0 * sin(x)) end - ff_harmonic! = DynamicalODEFunction(f1_harmonic!, f2_harmonic!; analytic = harmonic_analytic) + 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, Newmark(), dense_errors = true) + 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(du, u, p, t) - # return -u - 0.5 * du - # end - # function damped_oscillator!(ddu, du, u, p, t) - # @. ddu = -u - 0.5 * du - # return nothing - # end + 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) - OrdinaryDiffEq.SciMLBase.ArrayPartition( + 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)) @@ -42,40 +48,24 @@ end 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((ddu, v, u, p, t) -> ddu = -u - 0.5 * v, - (du, v, u, p, t) -> du = v, - analytic = (du0_u0, p, t) -> damped_oscillator_analytic ) + 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, Newmark(), dense_errors = true) + sim = test_convergence(dts, prob, NewmarkBeta(), dense_errors = true) @test sim.𝒪est[:l2]≈2 rtol=1e-1 -end - -# @testset "in-place vs. out-of-place" begin -# ode_i = SecondOrderODEProblem(damped_oscillator!, -# [0.0], [1.0], -# (0.0, 10.0)) -# ode_o = SecondOrderODEProblem(damped_oscillator, -# [0.0], [1.0], -# (0.0, 10.0)) -# @testset "NewmarkBeta" begin -# alg = NewmarkBeta() -# dt = 0.5 -# # fixed time step -# sol_i = solve(ode_i, alg, dt = dt) -# sol_o = solve(ode_o, alg, dt = dt) -# @test sol_i.t ≈ sol_o.t -# @test sol_i.u ≈ sol_o.u -# @test sol_i.destats.nf == sol_o.destats.nf -# @test sol_i.destats.nf2 == sol_o.destats.nf2 -# @test sol_i.destats.naccept == sol_o.destats.naccept -# @test 19 <= sol_i.destats.naccept <= 21 -# @test abs(sol_i.destats.nf - 4 * sol_i.destats.naccept) < 4 -# end -# end + # TODO + # function damped_oscillator(v, u, p, t) + # return -u - 0.5 * v + # end + # ... +end