Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Boundary conditions should always use solution object #208

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Boundary conditions should always use solution object #208

wants to merge 6 commits into from

Conversation

ErikQQY
Copy link
Member

@ErikQQY ErikQQY commented Aug 19, 2024

Fix #107
Fix #185

We should always pass a solution type to boundary conditions, and with this PR, the definition of boundary conditions for MIRK methods and Shooting Methods is consistent. I used linear interpolation to interpolate the solution, not sure it is optimal, but it is working fine. Now there is only one way to define our boundary condition:

using BoundaryValueDiffEq
tspan = (0.0, pi / 2)
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -9.81 * sin(θ)
end
function bc!(residual, u, p, t)
    residual[1] = u(pi/4)[1] + pi / 2
    residual[2] = u(pi/2)[1] - pi / 2
end
prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
sol = solve(prob, MIRK4(), dt = 0.05)

ErikQQY and others added 4 commits August 19, 2024 17:23
@ErikQQY
Boundary conditions should always use solution object
0454588 
Signed-off-by: ErikQQY <2283984853@qq.com>
@ErikQQY
Fix linear interpolation error
1ac7636 
Signed-off-by: ErikQQY <2283984853@qq.com>
@ErikQQY
Update benchmark scripts
ad85703
@ErikQQY
Fix indexing with MIRKSol
fd2be9b 
Signed-off-by: ErikQQY <2283984853@qq.com>
ChrisRackauckas
ChrisRackauckas reviewed Aug 19, 2024
View reviewed changes
src/solve/mirk.jl Outdated Show resolved Hide resolved
@ChrisRackauckas
Copy link
Member

used linear interpolation to interpolate the solution, not sure it is optimal, but it is working fine

It's not optimal. Each tableau has an associated interpolation we should add to the solution.

@ErikQQY
Change to MIRK interpolation
7b9a514 
Signed-off-by: ErikQQY <2283984853@qq.com>
@github-actions GitHub Actions
Copy link
Contributor

Benchmark Results

master 7b9a514... master/7b9a514a49c678...
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK2() 6.54 ± 0.2 ms 6.58 ± 0.16 ms 0.994
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK3() 2.35 ± 0.088 ms 2.41 ± 0.075 ms 0.976
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK4() 0.822 ± 0.032 ms 0.849 ± 0.026 ms 0.969
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK5() 0.869 ± 0.035 ms 0.897 ± 0.035 ms 0.968
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK6() 0.993 ± 0.033 ms 1.01 ± 0.03 ms 0.981
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false) 1.32 ± 0.1 ms 1.33 ± 0.095 ms 0.999
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true) 2.57 ± 0.21 ms 2.58 ± 0.16 ms 0.996
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.0354 ± 0.0048 s 0.0346 ± 0.0049 s 1.02
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.117 ± 0.0042 s 0.115 ± 0.0027 s 1.02
Simple Pendulum/IIP/Shooting(Tsit5()) 0.188 ± 0.0048 ms 0.186 ± 0.0062 ms 1.02
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK2() 0.0433 ± 0.0016 s 0.0426 ± 0.002 s 1.02
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK3() 12.4 ± 0.21 ms 12.2 ± 0.13 ms 1.01
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK4() 3.53 ± 0.077 ms 3.5 ± 0.063 ms 1.01
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK5() 3.64 ± 0.083 ms 3.57 ± 0.063 ms 1.02
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK6() 3.76 ± 0.091 ms 3.72 ± 0.082 ms 1.01
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false) 3.44 ± 0.88 ms 3.37 ± 0.81 ms 1.02
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true) 6.14 ± 1.3 ms 5.94 ± 1.3 ms 1.03
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.128 ± 0.0097 s 0.125 ± 0.012 s 1.02
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.4 ± 0.012 s 0.384 ± 0.0034 s 1.04
Simple Pendulum/OOP/Shooting(Tsit5()) 0.796 ± 0.036 ms 0.751 ± 0.045 ms 1.06
time_to_load 5.85 ± 0.083 s 5.91 ± 0.031 s 0.989

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

ChrisRackauckas
ChrisRackauckas reviewed Aug 19, 2024
View reviewed changes
src/solve/mirk.jl
recursive_flatten_twopoint!(resid, resids, cache.resid_size)
return nothing
end

@views function __mirk_loss(u, p, y, pt::StandardBVProblem, bc::BC, mesh, cache) where {BC}
y_ = recursive_unflatten!(y, u)
soly_ = MIRKSol(y_, mesh, cache)
resid_bc = eval_bc_residual(pt, bc, soly_, p, mesh)
soly_ = SciMLBase.build_solution(cache.prob, cache.alg, mesh, y_, interp = MIRKInterpolation(cache.mesh, y_, cache))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you should be able to build the solution type once and just keep updating its values?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, that's better

@ChrisRackauckas
Copy link
Member

This was done in a different PR right?

@ErikQQY
Copy link
Member Author

ErikQQY commented Sep 4, 2024

The changing of boundary condition definitions from res[1] = u[:, end] - 1.0 to res[1] = u(pi/2) - 1.0 will be done in this PR, though there are still some issues with building the solution type in the begining, this will be done very soon.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@ChrisRackauckas ChrisRackauckas ChrisRackauckas left review comments

Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

Indexing in boundary conditions README example should use interpolation instead of indexing
2 participants
@ErikQQY @ChrisRackauckas