Breaking Changes
The following git diff shows a comparison between the script in Appendix B of the original GridapTopOpt paper and the most up-to-date version of the package. Details regarding all breaking changes can be found in the next section.
using GridapTopOpt, Gridap
function main(write_dir)
# FE parameters
order = 1 # Finite element order
xmax,ymax = (1.0,1.0) # Domain size
dom = (0,xmax,0,ymax) # Bounding domain
el_size = (200,200) # Mesh partition size
prop_Γ_N = 0.2 # Γ_N size parameter
prop_Γ_D = 0.2 # Γ_D size parameter
f_Γ_N(x) = (x[1] ≈ xmax && # Γ_N indicator function
ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + eps())
f_Γ_D(x) = (x[1] ≈ 0.0 && # Γ_D indicator function
(x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps()))
# FD parameters
γ = 0.1 # HJ eqn time step coeff
γ_reinit = 0.5 # Reinit. eqn time step coeff
max_steps = floor(Int,order*minimum(el_size)/10)# Max steps for advection
tol = 1/(5order^2)/minimum(el_size) # Reinitialisation tolerance
# Problem parameters
κ = 1 # Diffusivity
g = 1 # Heat flow in
vf = 0.4 # Volume fraction constraint
lsf_func = initial_lsf(4,0.2) # Initial level set function
iter_mod = 10 # VTK Output modulo
path = "$write_dir/therm_comp_serial/" # Output path
mkpath(path) # Create path
# Model
model = CartesianDiscreteModel(dom,el_size);
update_labels!(1,model,f_Γ_D,"Gamma_D")
update_labels!(2,model,f_Γ_N,"Gamma_N")
# Triangulation and measures
Ω = Triangulation(model)
Γ_N = BoundaryTriangulation(model,tags="Gamma_N")
dΩ = Measure(Ω,2*order)
dΓ_N = Measure(Γ_N,2*order)
# Spaces
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"])
U = TrialFESpace(V,0.0)
V_φ = TestFESpace(model,reffe)
V_reg = TestFESpace(model,reffe;dirichlet_tags=["Gamma_N"])
U_reg = TrialFESpace(V_reg,0)
# Level set and interpolator
φh = interpolate(lsf_func,V_φ)
interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model)))
I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
# Weak formulation
- a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
- l(v,φ,dΩ,dΓ_N) = ∫(g*v)dΓ_N
+ a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+ l(v,φ) = ∫(g*v)dΓ_N
- state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+ state_map = AffineFEStateMap(a,l,U,V,V_φ)
# Objective and constraints
- J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
- dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+ J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+ dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
vol_D = sum(∫(1)dΩ)
- C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
- dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+ C1(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+ dC1(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
pcfs = PDEConstrainedFunctionals(J,[C1],state_map,
analytic_dJ=dJ,analytic_dC=[dC1])
# Velocity extension
α = 4max_steps*γ*maximum(get_el_Δ(model))
a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ
vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
# Finite difference scheme
scheme = FirstOrderStencil(length(el_size),Float64)
- ls_evo = HamiltonJacobiEvolution(scheme,model,V_φ,tol,max_steps)
+ evo = FiniteDifferenceEvolver(scheme,model,V_φ;max_steps)
+ reinit = FiniteDifferenceReinitialiser(scheme,model,V_φ;tol,γ_reinit)
+ ls_evo = LevelSetEvolution(evo,reinit)
# Optimiser
- optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,γ_reinit,
- verbose=true)
+ optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;γ,verbose=true)
# Solve
for (it,uh,φh) in optimiser
data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
write_history(path*"/history.txt",get_history(optimiser))
end
# Final structure
it = get_history(optimiser).niter; uh = get_state(pcfs)
writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,
"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
end
main("results/")Details of all breaking changes
This page will be updated in the event that a breaking change is introduced into the source code.
Updating from v0.2 to v0.3/v0.4
In v0.3 we added Zygote compatability for backwards differentiation in serial and parallel. In addition, in v0.4 we removed the requirement for φh to be passed to the constructors of StateMaps (excepting staggered maps). Below we list any breaking changes that will require changes to scripts implemented in v0.2:
StateMapsnow always differentiate onto the space of the primal variable. See #76 for details. This introduces a breaking API change asU_regandφhare removed from constructors ofStateMaps. E.g.,
becomesAffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
Backwards compatability has been added here to ensure that the old API produces a meaningful error.AffineFEStateMap(a,l,U,V,V_φ)- The way that we allocate vectors for distributed has been reworked. Now, we always create derivatives using
zero(<:FESpace), we then move this to the correct type of array when required. For example, if needing to interpolate onto a RHS vector (this doesn't have ghosts in distributed), we can use the functionality_interpolate_onto_rhs!. This change cleans up a lot of the constructors for<:StateMapandStateMapWithParam. Any code implemented that relies on the old approach for allocating vectors should be adjusted accordingly.
Updating from v0.1 to v0.2
In v0.2 we made several quality of life changes and enabled compatability with GridapEmbedded. Below we list any breaking changes that will require changes to scripts implemented in v0.1:
- Automatic differentiation capability has now been added to GridapDistributed. As a result, the
IntegrandWithMeasurestructure has been removed. In addition functionals previously required the measures to be passed as arguments, e.g.,
This is no longer required and the above should instead be written asJ(u,φ,dΩ,dΓ_N) = ∫(f(u,φ))dΩ + ∫(g(u,φ))dΓ_NJ(u,φ) = ∫(f(u,φ))dΩ + ∫(g(u,φ))dΓ_N