Newmark beta method #2187
Newmark beta method #2187
Conversation
src/alg_utils.jl
Outdated
| @@ -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. | |||
There was a problem hiding this comment.
Since the coefficients are part of the algorithm, can you calculate the order here?
There was a problem hiding this comment.
similar to FBDF, you can just set this to 1 I think
There was a problem hiding this comment.
Thanks for taking a look Hendrik! Yes, I could compute them here. I will fix such stuff after I get a pendulum example working. I just opened the PR to discuss the progress and possibly better design choices for the new class of solvers.
|
The last remaining problem is a bit tricky. For a given DynamicalODEFunction Newmark needs to solve for a "future accelleration" by solving the problem The first issue I run into is passing down the analytical Jacobian: using OrdinaryDiffEq
function f1_harmonic(dv, v, u, p, t)
dv .= -u
end
function f2_harmonic(du, v, u, p, t)
du .= v
end
harmonic_jac_f1(J, v, u, p, t) = J[1,1] = -1.0
ff_harmonic = DynamicalODEFunction(f1_harmonic, f2_harmonic; jac = harmonic_jac_f1)
DiffEqBase.has_jac(ff_harmonic) # false
DiffEqBase.has_jac(ff_harmonic.f1) # false
ff_f1 = ODEFunction{false}(ODEFunction{false}(f1_harmonic); jac=harmonic_jac_f1!)
ff_harmonic = DynamicalODEFunction(ff_f1, f2_harmonic; jac = harmonic_jac_f1)
DiffEqBase.has_jac(ff_harmonic) # false
DiffEqBase.has_jac(ff_harmonic.f1) # false The second problem is getting the AD on track. My current strategy is to pass a helper type |
|
I think I can call this now a first prototype. I am not very sure why it manages to converge to the correct solution, but I assume the residual is just wrong up to scaling. Another problem is that the Jacobian is per construction wrong. The derivative of However, this seems to break a few assumptions taken in the code. Hence I could need some feedback on how to proceed before putting time into a design which we do not want in the library. I already think that Furthermore I cannot figure out how to integrate the W-transformation here. Any pointers? |
src/integrators/integrator_utils.jl
Outdated
| @@ -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 | |||
There was a problem hiding this comment.
I'm having a very hard time figuring out whether this seems obviously fine, or an incredibly bad idea.
There was a problem hiding this comment.
Yea, as you pointed out on Slack we should switch to SecondOrderODEFunction IIUC.
There was a problem hiding this comment.
Okay the reason why I did this was that there is not SecondOrderODEFunction. Can we easily infer from the DynamicalODEFunction that the second function is only f2(du,u,v,p,t) = v?
There was a problem hiding this comment.
oh, that's weird. we have a SecondOrderODEProblem but no SecondOrderODEFunction.
src/nlsolve/newton.jl
Outdated
| @@ -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))) | |||
There was a problem hiding this comment.
why is this random line commented out?
There was a problem hiding this comment.
I could not get it to work with Newmark, because here we have a comparison between a scalar W_γdt and a pair γdt.
src/nlsolve/newton.jl
Outdated
| 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 |
There was a problem hiding this comment.
why isn't this compute_ustep!?
There was a problem hiding this comment.
I need to double check what exactly lead to this change, but I think we do not need this anymore.
src/nlsolve/type.jl
Outdated
| 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 |
There was a problem hiding this comment.
if you're going to do this, I think it's probably worth refactoring to use innertmp and outertmp for consistency (and then we can unify the handling by setting the appropriate one to 0 for the method being used.
There was a problem hiding this comment.
I think this should go into a different PR.
| integrator.f(cache.fsalfirst, integrator.uprev, integrator.p, integrator.t) | ||
| integrator.stats.nf += 1 | ||
| integrator.fsalfirst = cache.fsalfirst | ||
| integrator.fsallast = cache.fsalfirst |
There was a problem hiding this comment.
If you ask like this, then I think I am understanding something wrong here. From looking at other FSAL methods I thought that we can just put the same memory into both locations.
There was a problem hiding this comment.
never mind on this one then. I haven't read through other ones in a while, so this looked like a copy paste error, but if other methods do it also, then I think this is fine.
| upred_full = ArrayPartition( | ||
| duprev + dt*(1.0 - γ)*dduprev, | ||
| uprev + dt*dt*(0.5 - β)*dduprev + dt*duprev | ||
| ) |
There was a problem hiding this comment.
this seems like it will allocate a bunch, right?
There was a problem hiding this comment.
Indeed. As also pointed out below by you during review we might be able to use nlsolver.tmp to get rid of this.
| # 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 |
There was a problem hiding this comment.
I believe nlsolver.tmp isn't modified, but I'm not 100% sure.
| # 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ₙ₊₁)'': |
There was a problem hiding this comment.
I think we generally put the Δtₙ as part of γ rather than z. Is there a reason not to do so?
There was a problem hiding this comment.
Note that the
There was a problem hiding this comment.
Thanks for the feedback. I have updated the implementation to master and made the tests for linear problems (harmonic and damped oscillator) pass. However, I had to manually unroll the implementation. I left some inline comments here on Github on the code sections where I am stuck (and on which I will resume working after fixing the most pressing issues).
There was a problem hiding this comment.
I will strip down the imports to the used ones after we have fixed the remaining parts.
| @@ -0,0 +1,5 @@ | |||
| alg_extrapolates(alg::NewmarkBeta) = true | |||
There was a problem hiding this comment.
I will recover the extrapolation once we have ironed out the remaining parts.
| 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") |
There was a problem hiding this comment.
@oscardssmith I think I need some help here regarding the integration with OrdinaryDiffEqNonlinearSolve.jl .
A critical missing piece is an interface for the evaluation of the Jacobian in a suitable form (especially via ad). For now I have just bypassed this with a custom jacobian function, e.g. for the harmonic oscillator in the test
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
because many implicit second order ODE methods requires that J = Δtₙ²β ∂fᵤ + Δtₙγ ∂fᵥ . See above for a "full" derivation.
| function f2_harmonic!(du, v, u, p, t) | ||
| du .= v | ||
| end |
There was a problem hiding this comment.
What is the intended way to enforce this automatically? We do not have a SecondOrderODEFunction analogue to DynamicalODEFunction.
Prototype to resolve #411 .
Checklist
contributor guidelines, in particular the SciML Style Guide and
COLPRAC.