Newmark beta method #2187

termi-official · 2024-05-10T18:49:38Z

Prototype to resolve #411 .

Checklist

Appropriate tests were added
Any code changes were done in a way that does not break public API
All documentation related to code changes were updated
The new code follows the
contributor guidelines, in particular the SciML Style Guide and
COLPRAC.
Any new documentation only uses public API

ranocha · 2024-05-10T19:04:30Z

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.


Since the coefficients are part of the algorithm, can you calculate the order here?

similar to FBDF, you can just set this to 1 I think

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.

termi-official · 2024-05-15T17:24:09Z

The last remaining problem is a bit tricky. For a given DynamicalODEFunction Newmark needs to solve for a "future accelleration" by solving the problem $Δtₙ f₁(innertmp₂ + z β Δtₙ², innertmp₁ + z γ Δtₙ,t) = z$ for $z$ with given $innertmp, \alpha,\beta$.

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 ArrayPartitionNLSolveHelper around, which constructs an ArrayPartition when multiplied with a Vector. The idea is to provide an interface for inplace evaluation of f.f1 (DynamicalODEFunction) which is compatible with the evaluation form required by the internal nonlinear solver, such that we can call $f₁(innertmp + z \alpha,t)$, where $\alpha$ is the mentioned ArrayPartitionNLSolveHelper. Should we head down this road or is there possibly a better design?

termi-official · 2024-05-15T22:48:05Z

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 $f_1$ w.r.t. the input parameter is a rectangular matrix. For Newmark we assume that the input space is partitioned into 2 equally sized parts (velocity $v$ and displacement $u$) and that $f_2$ is just $du = v$. The first assumption tells us that the Jacobian is blocked with 2 blocks. For the inner nonlinear solve in Newmark we now obtain the system matrix for the Newton as $\partial_z f_1(u+\gamma_1 z,v+\gamma_2 z) = \gamma_1 J_{1,1}(u+\gamma_1 z,v+\gamma_2 z) + \gamma_2 J_{1,2}(u+\gamma_1 z,v+\gamma_2 z)$ where $J_{i,j}$ is the i,j-th block of the Jacobian of $f$.

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 ArrayPartitionNLSolveHelper is not great, but I cannot figure out a better design by myself either.

Furthermore I cannot figure out how to integrate the W-transformation here. Any pointers?

oscardssmith · 2024-05-22T12:53:33Z

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


I'm having a very hard time figuring out whether this seems obviously fine, or an incredibly bad idea.

Yea, as you pointed out on Slack we should switch to SecondOrderODEFunction IIUC.

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?

oh, that's weird. we have a SecondOrderODEProblem but no SecondOrderODEFunction.

oscardssmith · 2024-05-22T12:58:43Z

src/misc_utils.jl

+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
+
+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
+    ArrayPartitionNLSolveHelper(scalar*γ.γ₁, scalar*γ.γ₂)
+end
+
+Base.inv(γ::ArrayPartitionNLSolveHelper) = 1/(γ.γ₁+γ.γ₂)


This seems like a really weird algebra to me. Why does inv(::ArrayPartitionNLSolveHelper) give a scalar? also having the multiplication with Vector seems very sketchy. Probably should be AbstractVector? but It also feels somewhat weird that that method returns an ArrayPartition?...

Because I am not sure how to integrate a pair of weights with the W transformation. This is basically a hotfix to make the example work at all.

If you take another look at the perform step function you will see what I am trying to achive. Basically the nonlinear solver interface assumes that $\gamma$ is a scalar, but for Newmark methods we have one gamme for the first part of the solution (block 1 in the associated ArrayPartition) and one gamma for block 2).

oscardssmith · 2024-05-22T12:59:31Z

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)))


why is this random line commented out?

I could not get it to work with Newmark, because here we have a comparison between a scalar W_γdt and a pair γdt.

oscardssmith · 2024-05-22T13:02:17Z

src/nlsolve/newton.jl

+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


why isn't this compute_ustep!?

I need to double check what exactly lead to this change, but I think we do not need this anymore.

oscardssmith · 2024-05-22T13:05:01Z

src/nlsolve/type.jl

+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


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.

I think this should go into a different PR.

oscardssmith · 2024-05-22T13:06:52Z