From d1ff2ce7dcd759c0ae1710173c9752d1558e2ec5 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Fri, 5 Jul 2024 13:33:07 +1000 Subject: [PATCH 1/5] Correction to file names --- scripts/Benchmarks/generate_benchmark_scripts.jl | 2 +- scripts/Benchmarks/jobtemplate.sh | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/Benchmarks/generate_benchmark_scripts.jl b/scripts/Benchmarks/generate_benchmark_scripts.jl index 071177a..4051e35 100644 --- a/scripts/Benchmarks/generate_benchmark_scripts.jl +++ b/scripts/Benchmarks/generate_benchmark_scripts.jl @@ -96,7 +96,7 @@ phys_types = [ ]; # Template -template = read("$(ENV["SCRATCH"])/$dir_name/scripts/Benchmarks/jobtemplate_gadi.sh",String) +template = read("$(ENV["SCRATCH"])/$dir_name/scripts/Benchmarks/jobtemplate.sh",String) ## Generate Jobs jobs_by_phys = map(x->(x[1],generate_jobs(template,x[1],x[2],x[3])),phys_types); diff --git a/scripts/Benchmarks/jobtemplate.sh b/scripts/Benchmarks/jobtemplate.sh index 4bee92f..ed24050 100644 --- a/scripts/Benchmarks/jobtemplate.sh +++ b/scripts/Benchmarks/jobtemplate.sh @@ -1,7 +1,7 @@ #!/bin/bash #PBS -P np01 -#PBS -q normal +#PBS -q normal #PBS -N "{{:name}}" #PBS -l ncpus={{:ncpus}} #PBS -l mem={{:mem}}GB @@ -12,7 +12,7 @@ source $HOME/hpc-environments/gadi/load-intel.sh PROJECT_DIR=$SCRATCH/{{:dir_name}}/ mpiexec -n {{:ncpus}} julia --project=$PROJECT_DIR --check-bounds no -O3 \ - $PROJECT_DIR/scripts/Benchmarks/benchmark_gadi.jl \ + $PROJECT_DIR/scripts/Benchmarks/benchmark.jl \ {{:name}} \ {{{:write_dir}}} \ {{:type}} \ 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 2/5] 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 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 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 3/5] 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 e307851..3d436f2 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. From 0373fe725dfd885b66471cdb68700786ed05da05 Mon Sep 17 00:00:00 2001 From: Z J Wegert <60646897+zjwegert@users.noreply.github.com> Date: Mon, 12 Aug 2024 13:47:07 +1000 Subject: [PATCH 4/5] add logo --- README.md | 2 +- assets/logo.svg | 116 +++++++++++++++++++++++++++++++++++++++ docs/src/assets/logo.svg | 116 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 233 insertions(+), 1 deletion(-) create mode 100644 assets/logo.svg create mode 100644 docs/src/assets/logo.svg diff --git a/README.md b/README.md index 396fe2b..0056246 100755 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# GridapTopOpt +# drawing

GridapTopOpt [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://zjwegert.github.io/GridapTopOpt.jl/stable) [![Build Status](https://github.com/zjwegert/GridapTopOpt.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/zjwegert/GridapTopOpt.jl/actions) [![Codecov](https://codecov.io/gh/zjwegert/GridapTopOpt.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/zjwegert/GridapTopOpt.jl) diff --git a/assets/logo.svg b/assets/logo.svg new file mode 100644 index 0000000..eace389 --- /dev/null +++ b/assets/logo.svg @@ -0,0 +1,116 @@ + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg new file mode 100644 index 0000000..eace389 --- /dev/null +++ b/docs/src/assets/logo.svg @@ -0,0 +1,116 @@ + + + + + + + + + + + + + + + + + + From 965a3f12758513ccc1f3db3442733f70464eca2d Mon Sep 17 00:00:00 2001 From: Z J Wegert <60646897+zjwegert@users.noreply.github.com> Date: Tue, 13 Aug 2024 11:46:52 +1000 Subject: [PATCH 5/5] Update README.md Update to show dev docs --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 0056246..61d9caf 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # drawing

GridapTopOpt -[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://zjwegert.github.io/GridapTopOpt.jl/stable) [![Build Status](https://github.com/zjwegert/GridapTopOpt.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/zjwegert/GridapTopOpt.jl/actions) [![Codecov](https://codecov.io/gh/zjwegert/GridapTopOpt.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/zjwegert/GridapTopOpt.jl) +[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://zjwegert.github.io/GridapTopOpt.jl/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://zjwegert.github.io/GridapTopOpt.jl/dev) [![Build Status](https://github.com/zjwegert/GridapTopOpt.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/zjwegert/GridapTopOpt.jl/actions) [![Codecov](https://codecov.io/gh/zjwegert/GridapTopOpt.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/zjwegert/GridapTopOpt.jl) GridapTopOpt is computational toolbox for level set-based topology optimisation implemented in Julia and the [Gridap](https://github.com/gridap/Gridap.jl) package ecosystem. See the documentation and following publication for further details: