From 6528d1e2b6c148ecafde2fdaff3cf358f8c62bf5 Mon Sep 17 00:00:00 2001
From: zjwegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 26 Jul 2024 10:10:59 +1000
Subject: [PATCH 1/2] Improve HPM efficiency for analytic sens cases

---
 src/Optimisers/HilbertianProjection.jl | 108 ++++++++++++++++++-----
 test/seq/ThermalComplianceHPMTests.jl  | 114 +++++++++++++++++++++++++
 test/seq/runtests.jl                   |   1 +
 3 files changed, 201 insertions(+), 22 deletions(-)
 create mode 100644 test/seq/ThermalComplianceHPMTests.jl

diff --git a/src/Optimisers/HilbertianProjection.jl b/src/Optimisers/HilbertianProjection.jl
index 94d8033f..e307851c 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 00000000..87fce14c
--- /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 e89e16e1..595f3f10 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

From c327239c1cd3c22e90166bc7636a98ea023abd9e Mon Sep 17 00:00:00 2001
From: Z J Wegert <60646897+zjwegert@users.noreply.github.com>
Date: Fri, 26 Jul 2024 10:47:38 +1000
Subject: [PATCH 2/2] Typo in comment

---
 src/Optimisers/HilbertianProjection.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Optimisers/HilbertianProjection.jl b/src/Optimisers/HilbertianProjection.jl
index e307851c..3d436f2b 100644
--- a/src/Optimisers/HilbertianProjection.jl
+++ b/src/Optimisers/HilbertianProjection.jl
@@ -222,7 +222,7 @@ struct HilbertianProjection{A} <: Optimiser
 
   !!! note
       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
+      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.