diff --git a/Project.toml b/Project.toml index 2410fd74ca..365544eacb 100644 --- a/Project.toml +++ b/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" 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/Project.toml b/lib/OrdinaryDiffEqNewmark/Project.toml new file mode 100644 index 0000000000..a829d2589f --- /dev/null +++ b/lib/OrdinaryDiffEqNewmark/Project.toml @@ -0,0 +1,46 @@ +name = "OrdinaryDiffEqNewmark" +uuid = "d204908a-63b9-11ef-18f5-cfec7123a93b" +authors = ["Dennis Ogiermann "] +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"] diff --git a/lib/OrdinaryDiffEqNewmark/src/OrdinaryDiffEqNewmark.jl b/lib/OrdinaryDiffEqNewmark/src/OrdinaryDiffEqNewmark.jl new file mode 100644 index 0000000000..c773b4c86c --- /dev/null +++ b/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 \ No newline at end of file diff --git a/lib/OrdinaryDiffEqNewmark/src/alg_utils.jl b/lib/OrdinaryDiffEqNewmark/src/alg_utils.jl new file mode 100644 index 0000000000..a83c5e35b2 --- /dev/null +++ b/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 diff --git a/lib/OrdinaryDiffEqNewmark/src/algorithms.jl b/lib/OrdinaryDiffEqNewmark/src/algorithms.jl new file mode 100644 index 0000000000..9a22a70758 --- /dev/null +++ b/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 diff --git a/lib/OrdinaryDiffEqNewmark/src/newmark_caches.jl b/lib/OrdinaryDiffEqNewmark/src/newmark_caches.jl new file mode 100644 index 0000000000..d536ce7dda --- /dev/null +++ b/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 diff --git a/lib/OrdinaryDiffEqNewmark/src/newmark_nlsolve.jl b/lib/OrdinaryDiffEqNewmark/src/newmark_nlsolve.jl new file mode 100644 index 0000000000..e69de29bb2 diff --git a/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl b/lib/OrdinaryDiffEqNewmark/src/newmark_perform_step.jl new file mode 100644 index 0000000000..d01f736803 --- /dev/null +++ b/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 diff --git a/lib/OrdinaryDiffEqNewmark/test/runtests.jl b/lib/OrdinaryDiffEqNewmark/test/runtests.jl new file mode 100644 index 0000000000..4b90eb2815 --- /dev/null +++ b/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 diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/nlsolve.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/nlsolve.jl index 1633531bdc..235c17d053 100644 --- a/lib/OrdinaryDiffEqNonlinearSolve/src/nlsolve.jl +++ b/lib/OrdinaryDiffEqNonlinearSolve/src/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::NL, integrator::DiffEqBase.DEIntegrator, cache = nothing, repeat_step = false) where {NL <: AbstractNLSolver} diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl index 4c36666b6e..b6a334ca1f 100644 --- a/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl +++ b/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl @@ -72,7 +72,7 @@ mutable struct NLSolver{algType, iip, uType, gamType, tmpType, tType, C <: AbstractNLSolverCache, E} <: AbstractNLSolver{algType, iip} z::uType tmp::uType # DIRK and multistep methods only use tmp - tmp2::tmpType # for GLM if necessary + tmp2::tmpType # for GLM if neccssary ztmp::uType γ::gamType c::tType @@ -140,7 +140,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/lib/OrdinaryDiffEqRKN/test/nystrom_convergence_tests.jl b/lib/OrdinaryDiffEqRKN/test/nystrom_convergence_tests.jl index 107cf34235..f18ce6fc81 100644 --- a/lib/OrdinaryDiffEqRKN/test/nystrom_convergence_tests.jl +++ b/lib/OrdinaryDiffEqRKN/test/nystrom_convergence_tests.jl @@ -15,6 +15,8 @@ end ff_harmonic = DynamicalODEFunction(f1_harmonic, f2_harmonic; analytic = harmonic_analytic) prob = DynamicalODEProblem(ff_harmonic, v0, u0, (0.0, 5.0)) + + # Methods need BigFloat to test convergence rate dts = big"1.0" ./ big"2.0" .^ (5:-1:1) prob_big = DynamicalODEProblem(ff_harmonic, [big"1.0", big"1.0"], diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 503484ac33..82c52bd9bc 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -206,6 +206,9 @@ export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3B, EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43 +using OrdinaryDiffEqNewmark +export NewmarkBeta + import PrecompileTools import Preferences import DocStringExtensions