Added some serial tests + bug fix in ChainRules

scripts/MPI/3d_thermal_compliance_ALM.jl

@@ -29,7 +29,7 @@ function main(mesh_partition,distribute,el_size)
   dom = (0,xmax,0,ymax,0,zmax)
   γ = 0.1
   γ_reinit = 0.5
-  max_steps = floor(Int,order*minimum(el_size)/step)
+  max_steps = floor(Int,order*minimum(el_size)/10)
   tol = 1/(coef*order^2)/minimum(el_size)
   κ = 1
   η_coeff = 2

src/ChainRules.jl

@@ -65,13 +65,15 @@ Gridap.gradient(F::IntegrandWithMeasure,uh) = Gridap.gradient(F,[uh],1)
 Given an an `IntegrandWithMeasure` `F` and a vector of `FEFunctions` or `CellField` `uh` 
 (excluding measures) evaluate the Jacobian `F.F` with respect to `uh[K]`.
 """
-function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:FEFunction},K::Int)
+function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:Union{FEFunction,CellField}},K::Int)
   @check 0 < K <= length(uh)
   _f(uk) = F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...)
   return Gridap.jacobian(_f,uh[K])
 end
 
-function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector,K::Int)
+DistributedFields = Union{DistributedCellField,DistributedMultiFieldCellField}
+
+function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:DistributedFields},K::Int)
   @check 0 < K <= length(uh)
   local_fields = map(local_views,uh) |> to_parray_of_arrays
   local_measures = map(local_views,F.dΩ) |> to_parray_of_arrays

src/LevelSetTopOpt.jl

@@ -21,8 +21,8 @@ using Gridap: writevtk
 using GridapDistributed
 using GridapDistributed: DistributedDiscreteModel, DistributedTriangulation, 
   DistributedFESpace, DistributedDomainContribution, to_parray_of_arrays,
-  allocate_in_domain, DistributedCellField, DistributedMultiFieldFEBasis,
-  BlockPMatrix, BlockPVector, change_ghost
+  allocate_in_domain, DistributedCellField, DistributedMultiFieldCellField,
+  DistributedMultiFieldFEBasis, BlockPMatrix, BlockPVector, change_ghost
 
 using PartitionedArrays
 using PartitionedArrays: getany, tuple_of_arrays, matching_ghost_indices

test/seq/InverseHomogenisationALMTests.jl

@@ -0,0 +1,103 @@
+module InverseHomogenisationALMTests
+using Test
+
+using Gridap, LevelSetTopOpt
+
+"""
+  (Serial) Maximum bulk modulus inverse homogenisation with augmented Lagrangian method in 2D.
+
+  Optimisation problem:
+      Min J(Ω) = -κ(Ω)
+        Ω
+    s.t., Vol(Ω) = vf,
+          ⎡For unique εᴹᵢ, find uᵢ∈V=H¹ₚₑᵣ(Ω)ᵈ, 
+          ⎣∫ ∑ᵢ C ⊙ ε(uᵢ) ⊙ ε(vᵢ) dΩ = ∫ -∑ᵢ C ⊙ ε⁰ᵢ ⊙ ε(vᵢ) dΩ, ∀v∈V.
+""" 
+function main(;AD)
+  ## Parameters
+  order = 1
+  xmax,ymax=(1.0,1.0)
+  dom = (0,xmax,0,ymax)
+  el_size = (20,20)
+  γ = 0.05
+  γ_reinit = 0.5
+  max_steps = floor(Int,order*minimum(el_size)/10)
+  tol = 1/(5order^2)/minimum(el_size)
+  C = isotropic_elast_tensor(2,1.,0.3)
+  η_coeff = 2
+  α_coeff = 4max_steps*γ
+  vf = 0.5
+
+  ## FE Setup
+  model = CartesianDiscreteModel(dom,el_size,isperiodic=(true,true))
+  el_Δ = get_el_Δ(model)
+  f_Γ_D(x) = iszero(x)
+  update_labels!(1,model,f_Γ_D,"origin")
+
+  ## Triangulations and measures
+  Ω = Triangulation(model)
+  dΩ = Measure(Ω,2*order)
+  vol_D = sum(∫(1)dΩ)
+
+  ## Spaces
+  reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
+  reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe;dirichlet_tags=["origin"])
+  U = TrialFESpace(V,VectorValue(0.0,0.0))
+  V_reg = V_φ = TestFESpace(model,reffe_scalar)
+  U_reg = TrialFESpace(V_reg)
+
+  ## Create FE functions
+  lsf_fn = x->max(initial_lsf(2,0.4)(x),initial_lsf(2,0.4;b=VectorValue(0,0.5))(x))
+  φh = interpolate(lsf_fn,V_φ)
+
+  ## Interpolation and weak form
+  interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
+  I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
+
+  εᴹ = (TensorValue(1.,0.,0.,0.),     # ϵᵢⱼ⁽¹¹⁾≡ϵᵢⱼ⁽¹⁾
+        TensorValue(0.,0.,0.,1.),     # ϵᵢⱼ⁽²²⁾≡ϵᵢⱼ⁽²⁾
+        TensorValue(0.,1/2,1/2,0.))   # ϵᵢⱼ⁽¹²⁾≡ϵᵢⱼ⁽³⁾
+
+  a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ
+  l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
+
+  ## Optimisation functionals
+  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ))
+  dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ))
+  Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ,dΩ) = ∫(-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 = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ)
+  pcfs = if AD
+    PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol])
+  else
+    PDEConstrainedFunctionals(κ,[Vol],state_map)
+  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 = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
+    γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+
+  # Do a few iterations
+  vars, state = iterate(optimiser)
+  vars, state = iterate(optimiser,state)
+  true
+end
+
+# Test that these run successfully
+@test main(;AD=true)
+@test main(;AD=false)
+
+end # module

test/seq/InverterHPMTests.jl

@@ -0,0 +1,114 @@
+module InverterHPMTests
+using Test
+
+using Gridap, LevelSetTopOpt
+
+"""
+  (Serial) Inverter mechanism with Hilbertian projection method in 2D.
+
+  Optimisation problem:
+      Min J(Ω) = ηᵢₙ*∫ u⋅e₁ dΓᵢₙ/Vol(Γᵢₙ)
+        Ω
+    s.t., Vol(Ω) = vf,
+            C(Ω) = 0, 
+          ⎡u∈V=H¹(Ω;u(Γ_D)=0)ᵈ, 
+          ⎣∫ C ⊙ ε(u) ⊙ ε(v) dΩ + ∫ kₛv⋅u dΓₒᵤₜ = ∫ v⋅g dΓᵢₙ , ∀v∈V.
+        
+    where C(Ω) = ∫ -u⋅e₁-δₓ dΓₒᵤₜ/Vol(Γₒᵤₜ). We assume symmetry in the problem to aid
+     convergence.
+""" 
+function main()
+  ## Parameters
+  order = 1
+  dom = (0,1,0,0.5)
+  el_size = (20,20)
+  γ = 0.1
+  γ_reinit = 0.5
+  max_steps = floor(Int,order*minimum(el_size)/10)
+  tol = 1/(5order^2)/minimum(el_size)
+  C = isotropic_elast_tensor(2,1.0,0.3)
+  η_coeff = 2
+  α_coeff = 4max_steps*γ
+  vf = 0.4
+  δₓ = 0.2
+  ks = 0.1
+  g = VectorValue(0.5,0)
+
+  ## FE Setup
+  model = CartesianDiscreteModel(dom,el_size)
+  el_Δ = get_el_Δ(model)
+  f_Γ_in(x) = (x[1] ≈ 0.0) && (x[2] <= 0.03 + eps())
+  f_Γ_out(x) = (x[1] ≈ 1.0) && (x[2] <= 0.07 + eps())
+  f_Γ_D(x) = (x[1] ≈ 0.0) && (x[2] >= 0.4)
+  f_Γ_D2(x) = (x[2] ≈ 0.0)
+  update_labels!(1,model,f_Γ_in,"Gamma_in")
+  update_labels!(2,model,f_Γ_out,"Gamma_out")
+  update_labels!(3,model,f_Γ_D,"Gamma_D")
+  update_labels!(4,model,f_Γ_D2,"SymLine")
+
+  ## Triangulations and measures
+  Ω = Triangulation(model)
+  Γ_in = BoundaryTriangulation(model,tags="Gamma_in")
+  Γ_out = BoundaryTriangulation(model,tags="Gamma_out")
+  dΩ = Measure(Ω,2order)
+  dΓ_in = Measure(Γ_in,2order)
+  dΓ_out = Measure(Γ_out,2order)
+  vol_D = sum(∫(1)dΩ)
+  vol_Γ_in = sum(∫(1)dΓ_in)
+  vol_Γ_out = sum(∫(1)dΓ_out)
+
+  ## Spaces
+  reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
+  reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D","SymLine"],
+    dirichlet_masks=[(true,true),(false,true)])
+  U = TrialFESpace(V,[VectorValue(0.0,0.0),VectorValue(0.0,0.0)])
+  V_φ = TestFESpace(model,reffe_scalar)
+  V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_in","Gamma_out"])
+  U_reg = TrialFESpace(V_reg,[0,0])
+
+  ## Create FE functions
+  lsf_fn(x) = min(max(initial_lsf(6,0.2)(x),-sqrt((x[1]-1)^2+(x[2]-0.5)^2)+0.2),
+    sqrt((x[1])^2+(x[2]-0.5)^2)-0.1)
+  φh = interpolate(lsf_fn,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Γ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
+  l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in
+
+  ## Optimisation functionals
+  e₁ = VectorValue(1,0)
+  J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
+  Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
+
+  ## 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Γ_in,dΓ_out)
+  pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing])
+
+  ## 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,debug=true,constraint_names=[:Vol,:UΓ_out])
+
+  # Do a few iterations
+  vars, state = iterate(optimiser)
+  vars, state = iterate(optimiser,state)
+  true
+end
+
+# Test that these run successfully
+@test main()
+
+end # module

test/seq/NonlinearThermalComplianceALMTests.jl

@@ -0,0 +1,101 @@
+module NonlinearThermalComplianceALMTests
+using Test
+
+using Gridap, LevelSetTopOpt
+
+"""
+  (Serial) Minimum thermal compliance with augmented Lagrangian method in 2D with nonlinear diffusivity.
+
+  Optimisation problem:
+      Min J(Ω) = ∫ κ(u)*∇(u)⋅∇(u) dΩ
+        Ω
+    s.t., Vol(Ω) = vf,
+          ⎡u∈V=H¹(Ω;u(Γ_D)=0),
+          ⎣∫ κ(u)*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V.
+  
+  In this example κ(u) = κ0*(exp(ξ*u))
+""" 
+function main()
+  ## Parameters
+  order = 1
+  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/(5order^2)/minimum(el_size)
+  η_coeff = 2
+  α_coeff = 4max_steps*γ
+  vf = 0.4
+
+  ## 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.ρ
+
+  κ0 = 1
+  ξ = -1
+  κ(u) = κ0*(exp(ξ*u))
+  res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
+
+  ## Optimisation functionals
+  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))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 = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
+
+  ## 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 = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
+    γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+
+  # Do a few iterations
+  vars, state = iterate(optimiser)
+  vars, state = iterate(optimiser,state)
+  true
+end
+
+# Test that these run successfully
+@test main()
+
+end # module