diff --git a/src/Optimisers/HilbertianProjection.jl b/src/Optimisers/HilbertianProjection.jl index 94d8033..e307851 100644 --- a/src/Optimisers/HilbertianProjection.jl +++ b/src/Optimisers/HilbertianProjection.jl @@ -123,6 +123,13 @@ function compute_α(C,dC_orthog,dC,normsq,K,P,λ,α_min,α_max) 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 @@ -141,7 +148,7 @@ A Hilbertian projection method as described by Wegert et al., 2023 - `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 @@ -205,8 +212,8 @@ struct HilbertianProjection <: Optimiser - `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. @@ -214,9 +221,12 @@ struct HilbertianProjection <: Optimiser 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, especically + 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}, @@ -226,7 +236,7 @@ struct HilbertianProjection <: Optimiser 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, @@ -240,12 +250,15 @@ struct HilbertianProjection <: Optimiser 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 @@ -339,9 +352,26 @@ function Base.iterate(m::HilbertianProjection,state) 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 @@ -375,18 +405,52 @@ function Base.iterate(m::HilbertianProjection,state) ## 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, γ + 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) + end + end + + return J, C, dJ, dC, γ end function debug_print(orthog,dV,dC,dC_orthog,K,P,nullity,debug_code,debug) diff --git a/test/seq/ThermalComplianceHPMTests.jl b/test/seq/ThermalComplianceHPMTests.jl new file mode 100644 index 0000000..87fce14 --- /dev/null +++ b/test/seq/ThermalComplianceHPMTests.jl @@ -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 \ No newline at end of file diff --git a/test/seq/runtests.jl b/test/seq/runtests.jl index e89e16e..595f3f1 100644 --- a/test/seq/runtests.jl +++ b/test/seq/runtests.jl @@ -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