diff --git a/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl b/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl index cc671e2..6a9e2fd 100644 --- a/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl +++ b/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl @@ -1,7 +1,8 @@ module ODEProblemLibrary using DiffEqBase: DiffEqBase, ODEFunction, ODEProblem -using LinearAlgebra: LinearAlgebra, Diagonal, I, SymTridiagonal, diagind, ldiv!, mul!, rmul! +using LinearAlgebra: LinearAlgebra, Diagonal, I, SymTridiagonal, diagind, ldiv!, lmul!, mul!, + rmul! using Random: Random Random.seed!(100) diff --git a/lib/ODEProblemLibrary/src/filament_prob.jl b/lib/ODEProblemLibrary/src/filament_prob.jl index 04bd2fd..60bb243 100644 --- a/lib/ODEProblemLibrary/src/filament_prob.jl +++ b/lib/ODEProblemLibrary/src/filament_prob.jl @@ -227,13 +227,13 @@ function projection!(f::FilamentCache) mul!(J_JT, J, J') LDLt_inplace!(J_JT_LDLT, J_JT) ldiv!(P0, J_JT_LDLT, J) - mul!(P', P0, J) + mul!(P, P0', J) subtract_from_identity!(P) return nothing end function subtract_from_identity!(A) - rmul!(-1, A) + lmul!(-1, A) @inbounds for i in 1:size(A, 1) A[i, i] += 1 end diff --git a/lib/ODEProblemLibrary/test/filament.jl b/lib/ODEProblemLibrary/test/filament.jl new file mode 100644 index 0000000..0801021 --- /dev/null +++ b/lib/ODEProblemLibrary/test/filament.jl @@ -0,0 +1,34 @@ +using ODEProblemLibrary, LinearAlgebra, Test + +# The right hand side is dr = P(r) * (A * r + F(t)), where P projects onto the +# tangent space of the inextensibility constraints, so every squared segment +# length must have zero rate of change for any state and any time. +function max_segment_rate(prob, N, u, t) + du = similar(u) + prob.f(du, u, prob.p, t) + R, dR = reshape(u, 3, N + 1), reshape(du, 3, N + 1) + rate = maximum( + abs(2 * dot(R[:, i + 1] - R[:, i], dR[:, i + 1] - dR[:, i])) for i in 1:N + ) + return rate / max(1.0, maximum(abs, du)), du +end + +@testset "prob_ode_filament" begin + prob = prob_ode_filament + N = length(prob.u0) รท 3 - 1 + @test length(prob.u0) == 3 * (N + 1) + + for (u, t) in ( + (copy(prob.u0), 0.0), + (copy(prob.u0), 0.013), + (prob.u0 .+ 0.01 .* sinpi.(range(0, 2, length = length(prob.u0))), 0.013), + ) + rate, du = max_segment_rate(prob, N, u, t) + @test all(isfinite, du) + @test rate < 1.0e-10 + end + + J = zeros(length(prob.u0), length(prob.u0)) + prob.f.jac(J, prob.u0, prob.p, 0.0) + @test all(isfinite, J) +end diff --git a/lib/ODEProblemLibrary/test/runtests.jl b/lib/ODEProblemLibrary/test/runtests.jl index 8b9c556..e28caec 100644 --- a/lib/ODEProblemLibrary/test/runtests.jl +++ b/lib/ODEProblemLibrary/test/runtests.jl @@ -15,6 +15,7 @@ if TEST_GROUP == "Core" || TEST_GROUP == "All" @time @testset "Load Tests" begin @test ODEProblemLibrary isa Module end + @time @safetestset "Filament" include("filament.jl") end # Quality assurance: no undefined exports, stale dependencies, etc.