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.

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

Improve HPM efficiency #68

Merged
merged 2 commits into from
Aug 12, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 86 additions & 22 deletions src/Optimisers/HilbertianProjection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,13 @@
return α, ∑α², debug_flag
end

# Enable optimisations when not using auto diff. Without autodiff, we can compute
# sensitivities with a small cost and avoid recomputing the whole forward problem
# after the line search. When using AD we solve an adjoint problem for each constraint
# which is costly if done at each iteration of the line search.
struct WithAutoDiff end
struct NoAutoDiff end

"""
struct HilbertianProjection <: Optimiser

Expand All @@ -141,7 +148,7 @@
- `has_oscillations::Function`: A function to check for oscillations.
- `params::NamedTuple`: Optimisation parameters.
"""
struct HilbertianProjection <: Optimiser
struct HilbertianProjection{A} <: Optimiser
problem :: PDEConstrainedFunctionals
ls_evolver :: LevelSetEvolution
vel_ext :: VelocityExtension
Expand Down Expand Up @@ -205,18 +212,21 @@
- `ls_max_iters = 10`: Maximum number of line search iterations.
- `ls_δ_inc = 1.1`: Increase multiplier for `γ` on acceptance.
- `ls_δ_dec = 0.7`: Decrease multiplier for `γ` on rejection.
- `ls_ξ = 0.0025`: Line search tolerance for objective reduction.
- `ls_ξ_reduce_coef = 0.1`: Coeffient on `ls_ξ` if constraints within tolerance (see below).
- `ls_ξ = 1.0`: Line search tolerance for objective reduction.
- `ls_ξ_reduce_coef = 0.0025`: Coeffient on `ls_ξ` if constraints within tolerance (see below).
- `ls_ξ_reduce_abs_tol = 0.01`: Tolerance on constraints to reduce `ls_ξ` via `ls_ξ_reduce_coef`.
- `ls_γ_min = 0.001`: Minimum coeffient on the time step size for solving the HJ evolution equation.
- `ls_γ_max = 0.1`: Maximum coeffient on the time step size for solving the HJ evolution equation.

A more concervative evolution of the boundary can be achieved by decreasing `ls_γ_max`.

!!! note
For some problems (e.g., inverter mechanism), we have observed that a simple oscillation
detection algorithm leads to better convergence compared to the line search. By default
disabling the line search via `ls_enabled = false` will enable oscillation detection.
The line search has been adjusted so that it is only enforced once the constraints
are within a set tolerance. This generally leads to better optimisation histories, especially
for problems where constraints are far from saturation and the objective must decrease to
improve the constraints.

This can be set to always be enfored by taking `ls_ξ = 0.0025` and `ls_ξ_reduce_coef = 0.1`.
"""
function HilbertianProjection(
problem :: PDEConstrainedFunctionals{N},
Expand All @@ -226,7 +236,7 @@
orthog = HPModifiedGramSchmidt(),
λ=0.5, α_min=0.1, α_max=1.0, γ=0.1, γ_reinit=0.5, reinit_mod = 1,
ls_enabled = true, ls_max_iters = 10, ls_δ_inc = 1.1, ls_δ_dec = 0.7,
ls_ξ = 0.0025, ls_ξ_reduce_coef = 0.1, ls_ξ_reduce_abs_tol = 0.01,
ls_ξ = 1, ls_ξ_reduce_coef = 0.0025, ls_ξ_reduce_abs_tol = 0.01,
ls_γ_min = 0.001, ls_γ_max = 0.1,
maxiter = 1000, verbose=false, constraint_names = map(i -> Symbol("C_$i"),1:N),
converged::Function = default_hp_converged, debug = false,
Expand All @@ -240,12 +250,15 @@
al_keys = [:J,constraint_names...,:γ]
al_bundles = Dict(:C => constraint_names)
history = OptimiserHistory(Float64,al_keys,al_bundles,maxiter,verbose)

projector = HilbertianProjectionMap(N,orthog,vel_ext;λ,α_min,α_max,debug)

# Optimisation when not using AD
A = ((typeof(problem.analytic_dJ) <: Function) &&
all(@. typeof(problem.analytic_dC) <: Function)) ? NoAutoDiff : WithAutoDiff;

params = (;debug,γ,γ_reinit,reinit_mod,ls_enabled,ls_max_iters,ls_δ_inc,ls_δ_dec,ls_ξ,
ls_ξ_reduce_coef,ls_ξ_reduce_abs_tol,ls_γ_min,ls_γ_max,os_γ_mult)
new(problem,ls_evolver,vel_ext,projector,history,converged,
has_oscillations,params,φ0)
new{A}(problem,ls_evolver,vel_ext,projector,history,converged,has_oscillations,params,φ0)
end
end

Expand Down Expand Up @@ -339,9 +352,26 @@
U_reg = get_deriv_space(m.problem.state_map)
V_φ = get_aux_space(m.problem.state_map)
interpolate!(FEFunction(U_reg,θ),vel,V_φ)
J, C, dJ, dC, γ = _linesearch!(m,state)

## Hilbertian extension-regularisation
project!(m.vel_ext,dJ)
project!(m.vel_ext,dC)
θ = update_descent_direction!(m.projector,dJ,C,dC,m.vel_ext.K)

## Update history and build state
push!(history,(J,C...,γ))
uh = get_state(m.problem)
state = (;it=it+1, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, γ, os_it)
vars = params.debug ? (it,uh,φh,state) : (it,uh,φh)
return vars, state
end

function _linesearch!(m::HilbertianProjection{WithAutoDiff},state)
it, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, γ, os_it = state

ls_enabled = params.ls_enabled
reinit_mod = params.reinit_mod
params = m.params; history = m.history
ls_enabled = params.ls_enabled; reinit_mod = params.reinit_mod
ls_max_iters,δ_inc,δ_dec = params.ls_max_iters,params.ls_δ_inc,params.ls_δ_dec
ξ, ξ_reduce, ξ_reduce_tol = params.ls_ξ, params.ls_ξ_reduce_coef, params.ls_ξ_reduce_abs_tol
γ_min, γ_max = params.ls_γ_min,params.ls_γ_max
Expand Down Expand Up @@ -375,18 +405,52 @@

## Calculate objective, constraints, and shape derivatives after line search
J, C, dJ, dC = Gridap.evaluate!(m.problem,φh)
uh = get_state(m.problem)

## Hilbertian extension-regularisation
project!(m.vel_ext,dJ)
project!(m.vel_ext,dC)
θ = update_descent_direction!(m.projector,dJ,C,dC,m.vel_ext.K)
return J, C, dJ, dC, γ
end

## Update history and build state
push!(history,(J,C...,γ))
state = (;it=it+1, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, γ, os_it)
vars = params.debug ? (it,uh,φh,state) : (it,uh,φh)
return vars, state
function _linesearch!(m::HilbertianProjection{NoAutoDiff},state)
it, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, γ, os_it = state

params = m.params; history = m.history
ls_enabled = params.ls_enabled; reinit_mod = params.reinit_mod
ls_max_iters,δ_inc,δ_dec = params.ls_max_iters,params.ls_δ_inc,params.ls_δ_dec
ξ, ξ_reduce, ξ_reduce_tol = params.ls_ξ, params.ls_ξ_reduce_coef, params.ls_ξ_reduce_abs_tol
γ_min, γ_max = params.ls_γ_min,params.ls_γ_max

ls_it = 0; done = false
φ = get_free_dof_values(φh); copy!(φ_tmp,φ)
while !done && (ls_it <= ls_max_iters)
# Advect & Reinitialise
evolve!(m.ls_evolver,φ,vel,γ)
iszero(it % reinit_mod) && reinit!(m.ls_evolver,φ,params.γ_reinit)

# Check enabled
if ~ls_enabled
J, C, dJ, dC = Gridap.evaluate!(m.problem,φh)
return J, C, dJ, dC, γ

Check warning on line 431 in src/Optimisers/HilbertianProjection.jl

View check run for this annotation

Codecov / codecov/patch

src/Optimisers/HilbertianProjection.jl#L430-L431

Added lines #L430 - L431 were not covered by tests
end

# Calcuate new objective and constraints
J_interm, C_interm, dJ_interm, dC_interm = Gridap.evaluate!(m.problem,φh)

# Reduce line search parameter if constraints close to saturation
_ξ = all(Ci -> abs(Ci) < ξ_reduce_tol, C_interm) ? ξ*ξ_reduce : ξ

# Accept/reject
if (J_interm < J + _ξ*abs(J)) || (γ <= γ_min)
γ = min(δ_inc*γ, γ_max)
done = true
print_msg(history," Accepted iteration with γ = $(γ) \n";color=:yellow)
J, C, dJ, dC = J_interm, C_interm, dJ_interm, dC_interm
else
γ = max(δ_dec*γ, γ_min)
copy!(φ,φ_tmp)
print_msg(history," Reject iteration with γ = $(γ) \n";color=:red)

Check warning on line 449 in src/Optimisers/HilbertianProjection.jl

View check run for this annotation

Codecov / codecov/patch

src/Optimisers/HilbertianProjection.jl#L447-L449

Added lines #L447 - L449 were not covered by tests
end
end

return J, C, dJ, dC, γ
end

function debug_print(orthog,dV,dC,dC_orthog,K,P,nullity,debug_code,debug)
Expand Down
114 changes: 114 additions & 0 deletions test/seq/ThermalComplianceHPMTests.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
module ThermalComplianceHPMTests
using Test

using Gridap, GridapTopOpt
using GridapTopOpt: WithAutoDiff, NoAutoDiff
"""
(Serial) Minimum thermal compliance with augmented Lagrangian method in 2D.

Optimisation problem:
Min J(Ω) = ∫ κ*∇(u)⋅∇(u) dΩ
Ω
s.t., Vol(Ω) = vf,
⎡u∈V=H¹(Ω;u(Γ_D)=0),
⎣∫ κ*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V.
"""
function main(;order,AD_case)
## Parameters
xmax = ymax = 1.0
prop_Γ_N = 0.2
prop_Γ_D = 0.2
dom = (0,xmax,0,ymax)
el_size = (20,20)
γ = 0.1
γ_reinit = 0.5
max_steps = floor(Int,order*minimum(el_size)/10)
tol = 1/(5*order^2)/minimum(el_size)
κ = 1
vf = 0.4
η_coeff = 2
α_coeff = 4max_steps*γ

## FE Setup
model = CartesianDiscreteModel(dom,el_size);
el_Δ = get_el_Δ(model)
f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= ymax*prop_Γ_D + eps() ||
x[2] >= ymax-ymax*prop_Γ_D - eps()))
f_Γ_N(x) = (x[1] ≈ xmax && ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <=
ymax/2+ymax*prop_Γ_N/2 + eps())
update_labels!(1,model,f_Γ_D,"Gamma_D")
update_labels!(2,model,f_Γ_N,"Gamma_N")

## Triangulations and measures
Ω = Triangulation(model)
Γ_N = BoundaryTriangulation(model,tags="Gamma_N")
dΩ = Measure(Ω,2*order)
dΓ_N = Measure(Γ_N,2*order)
vol_D = sum(∫(1)dΩ)

## Spaces
reffe_scalar = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_D"])
U = TrialFESpace(V,0.0)
V_φ = TestFESpace(model,reffe_scalar)
V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"])
U_reg = TrialFESpace(V_reg,0)

## Create FE functions
φh = interpolate(initial_lsf(4,0.2),V_φ)

## Interpolation and weak form
interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ

a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N

## Optimisation functionals
J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ

## Finite difference solver and level set function
ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)

## Setup solver and FE operators
state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
pcfs = if AD_case == :no_ad
PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
elseif AD_case == :with_ad
PDEConstrainedFunctionals(J,[Vol],state_map)
elseif AD_case == :partial_ad1
PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ)
elseif AD_case == :partial_ad2
PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
else
@error "AD case not defined"
end

## Hilbertian extension-regularisation problems
α = α_coeff*maximum(el_Δ)
a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ;
vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)

## Optimiser
optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh;
γ,γ_reinit,verbose=true,constraint_names=[:Vol])

AD_case ∈ (:with_ad,:partial_ad1,:partial_ad2) && @test typeof(optimiser) <: HilbertianProjection{WithAutoDiff}
AD_case ∈ (:no_ad,) && @test typeof(optimiser) <: HilbertianProjection{NoAutoDiff}

# Do a few iterations
vars, state = iterate(optimiser)
vars, state = iterate(optimiser,state)
true
end

# Test that these run successfully
@test main(;order=1,AD_case=:with_ad)
@test main(;order=1,AD_case=:partial_ad1)
@test main(;order=1,AD_case=:partial_ad2)
@test main(;order=1,AD_case=:no_ad)

end # module
1 change: 1 addition & 0 deletions test/seq/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ module LSTOSequentialTests
using Test

@time @testset "Thermal Compliance - ALM" begin include("ThermalComplianceALMTests.jl") end
@time @testset "Thermal Compliance - HPM" begin include("ThermalComplianceHPMTests.jl") end
@time @testset "Nonlinear Thermal Compliance - ALM" begin include("NonlinearThermalComplianceALMTests.jl") end
@time @testset "Nonlinear Neohook with Jacobian - ALM" begin include("NeohookAnalyticJacALMTests.jl") end
@time @testset "Inverse Homogenisation - ALM" begin include("InverseHomogenisationALMTests.jl") end
Expand Down
Loading