Skip to content

Commit

Permalink
no allocs in state equation for trapeze (#320)
Browse files Browse the repository at this point in the history
  • Loading branch information
PierreMartinon authored Nov 15, 2024
1 parent 3dbd894 commit 34a3d2e
Show file tree
Hide file tree
Showing 2 changed files with 6 additions and 4 deletions.
4 changes: 2 additions & 2 deletions benchmark/prof.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ function local_mayer(obj, x0, xf, v)
end

function init(;grid_size, disc_method)
prob = goddard_all()
#prob = goddard()
#prob = goddard_all()
prob = goddard()
#prob = double_integrator_mintf()
ocp = prob[:ocp]
docp = CTDirect.DOCP(ocp, grid_size=grid_size, time_grid=CTDirect.__time_grid(), disc_method=string(disc_method))
Expand Down
6 changes: 4 additions & 2 deletions src/trapeze.jl
Original file line number Diff line number Diff line change
Expand Up @@ -116,8 +116,10 @@ function setConstraintBlock!(docp::DOCP{Trapeze}, c, xu, v, time_grid, i, work)
offset_dyn_i = (i-1)*docp.dim_NLP_x
offset_dyn_ip1 = i*docp.dim_NLP_x

# trapeze rule
@. c[offset+1:offset+docp.dim_OCP_x] = xip1 - (xi + 0.5 * (tip1 - ti) * (work[offset_dyn_i+1:offset_dyn_i+docp.dim_OCP_x] + work[offset_dyn_ip1+1:offset_dyn_ip1+docp.dim_OCP_x]))
# trapeze rule (no allocations ^^)
halfstep = 0.5 * (tip1 - ti)
@views @. c[offset+1:offset+docp.dim_OCP_x] = xip1 - (xi + halfstep * (work[offset_dyn_i+1:offset_dyn_i+docp.dim_OCP_x] + work[offset_dyn_ip1+1:offset_dyn_ip1+docp.dim_OCP_x]))

if docp.is_lagrange
c[offset+docp.dim_NLP_x] = get_lagrange_state_at_time_step(xu, docp, i+1) - (get_lagrange_state_at_time_step(xu, docp, i) + 0.5 * (tip1 - ti) * (work[offset_dyn_i+docp.dim_NLP_x] + work[offset_dyn_ip1+docp.dim_NLP_x]))
end
Expand Down

0 comments on commit 34a3d2e

Please sign in to comment.