added some testing

scripts/Embedded/Examples/2d_thermal_with_ad.jl

@@ -0,0 +1,95 @@
+using Pkg; Pkg.activate()
+
+using Gridap,GridapEmbedded,GridapTopOpt
+
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
+import Gridap.Geometry: get_node_coordinates, collect1d
+include("../embedded_measures.jl")
+
+function main()
+  path="./results/UnfittedFEM_thermal_compliance_ALM_AD/"
+  n = 100
+  order = 1
+  γ = 0.1
+  γ_reinit = 0.5
+  max_steps = floor(Int,order*minimum(n)/10)
+  tol = 1/(5*order^2)/minimum(n)
+  vf = 0.4
+  α_coeff = 4max_steps*γ
+  iter_mod = 1
+
+  model = CartesianDiscreteModel((0,1,0,1),(n,n));
+  el_Δ = get_el_Δ(model)
+  f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= 0.2 + eps() || x[2] >= 0.8 - eps()))
+  f_Γ_N(x) = (x[1] ≈ 1 && 0.4 - eps() <= x[2] <= 0.6 + 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)
+
+  φh = interpolate(x->-cos(4π*x[1])*cos(4*pi*x[2])/4-0.2/4,V_φ)
+  embedded_meas = EmbeddedMeasureCache(φh,V_φ)
+  update_meas(args...) = update_embedded_measures!(embedded_meas,args...)
+  get_meas(args...) = get_embedded_measures(embedded_meas,args...)
+
+  ## Weak form
+  a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-3*∇(u)⋅∇(v))dΩ2
+  l(v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(v)dΓ_N
+
+  ## Optimisation functionals
+  J(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u))dΩ1
+  Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf/vol_D)dΩ;
+
+  ## IntegrandWithEmbeddedMeasure
+  a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas)
+  l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),update_meas)
+  J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),update_meas)
+  Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),update_meas)
+
+  ## Evolution Method
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
+
+  ## Setup solver and FE operators
+  state_map = AffineFEStateMap(a_iem,l_iem,U,V,V_φ,U_reg,φh,(dΩ,dΓ_N))
+  pcfs = PDEConstrainedFunctionals(J_iem,[Vol_iem],state_map)
+
+  ## 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
+  rm(path,force=true,recursive=true)
+  mkpath(path)
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;reinit_mod=5,
+    γ,γ_reinit,verbose=true,constraint_names=[:Vol])
+  for (it,uh,φh) in optimiser
+    _Ω1,_,_Γ = get_embedded_triangulations(embedded_meas,φh)
+    if iszero(it % iter_mod)
+      writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+      writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
+    end
+    write_history(path*"/history.txt",optimiser.history)
+  end
+  it = get_history(optimiser).niter; uh = get_state(pcfs)
+  _Ω1,_,_Γ = get_embedded_triangulations(embedded_meas,φh)
+  writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+  writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
+end
+
+main()

scripts/Embedded/Examples/2d_thermal_with_ad_MPI.jl

@@ -0,0 +1,101 @@
+using Gridap,GridapTopOpt
+
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
+import Gridap.Geometry: get_node_coordinates, collect1d
+using GridapDistributed, GridapPETSc, PartitionedArrays
+
+include("../embedded_measures.jl")
+
+function main(parts,distribute)
+  ranks = distribute(LinearIndices((prod(parts),)))
+
+  path="./results/UnfittedFEM_thermal_compliance_ALM_MPI_AD/"
+  n = 64
+  order = 2
+  γ = 0.1
+  γ_reinit = 0.5
+  max_steps = floor(Int,order*minimum(n)/10)
+  tol = 1/(5*order^2)/minimum(n)
+  κ = 1
+  vf = 0.4
+  α_coeff = 4max_steps*γ
+  iter_mod = 1
+  i_am_main(ranks) && rm(path,force=true,recursive=true)
+  i_am_main(ranks) && mkpath(path)
+
+  model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(n,n));
+  el_Δ = get_el_Δ(model)
+  f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= 0.2 + eps() || x[2] >= 0.8 - eps()))
+  f_Γ_N(x) = (x[1] ≈ 1 && 0.4 - eps() <= x[2] <= 0.6 + 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)
+
+  φh = interpolate(initial_lsf(4,0.2),V_φ)
+  embedded_meas = EmbeddedMeasureCache(φh,V_φ)
+  update_meas(args...) = update_embedded_measures!(embedded_meas,args...)
+  get_meas(args...) = get_embedded_measures(embedded_meas,args...)
+
+  ## Weak form
+  a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-3*∇(u)⋅∇(v))dΩ2
+  l(v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(v)dΓ_N
+
+  ## Optimisation functionals
+  J(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u))dΩ1
+  Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf/vol_D)dΩ;
+
+  ## IntegrandWithEmbeddedMeasure
+  a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas)
+  l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),update_meas)
+  J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),update_meas)
+  Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),update_meas)
+
+  ## Evolution Method
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
+
+  ## Setup solver and FE operators
+  state_map = AffineFEStateMap(a_iem,l_iem,U,V,V_φ,U_reg,φh,(dΩ,dΓ_N))
+  pcfs = PDEConstrainedFunctionals(J_iem,[Vol_iem],state_map)
+
+  ## 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_mod=5,
+    γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
+  for (it,uh,φh) in optimiser
+    _Ω1,_,_Γ = get_embedded_triangulations(embedded_meas,φh)
+    if iszero(it % iter_mod)
+      writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+      writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
+    end
+    write_history(path*"/history.txt",optimiser.history;ranks)
+  end
+  it = get_history(optimiser).niter; uh = get_state(pcfs)
+  _Ω1,_,_Γ = get_embedded_triangulations(embedded_meas,φh)
+  writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+  writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
+end
+
+with_mpi() do distribute
+  parts = (3,3);
+  main(parts,distribute)
+end

scripts/Embedded/MWEs/gridap_embedded_grad_testing_dist_MWE.jl

@@ -11,7 +11,7 @@ import Gridap.Geometry: get_node_coordinates, collect1d
 
 include("../embedded_measures.jl")
 
-parts = (2,2);
+parts = (3,3);
 ranks = with_debug() do distribute
   distribute(LinearIndices((prod(parts),)))
 end
@@ -25,13 +25,31 @@ V = TestFESpace(model,reffe_scalar)
 U = TrialFESpace(V)
 V_φ = TestFESpace(model,reffe_scalar)
 
-φh = interpolate(x->-cos(4π*x[1])*cos(4*pi*x[2])/4-0.2/4,V_φ)
+# φh = interpolate(x->-cos(4π*x[1])*cos(4*pi*x[2])/4-0.2/4,V_φ)
+φh = interpolate(x->-sqrt((x[1]-1/2)^2+(x[2]-1/2)^2)+0.2,V_φ)
 embedded_meas = EmbeddedMeasureCache(φh,V_φ)
 update_meas(args...) = update_embedded_measures!(embedded_meas,args...)
 get_meas(args...) = get_embedded_measures(embedded_meas,args...)
+update_embedded_measures!(embedded_meas,φh)
 
 _j(φ,dΩ,dΩ1,dΩ2,dΓ) = ∫(φ)dΩ1 + ∫(φ)dΩ2 + ∫(φ)dΓ
 
 j_iem = IntegrandWithEmbeddedMeasure(_j,(dΩ,),update_meas)
 j_iem(φh)
-gradient(j_iem,φh)
+gradient(j_iem,φh)
+
+# Possible fix - doesn't work though!
+# function alt_get_geo_params(ϕh::CellField,Vbg,gids)
+#   Ωbg = get_triangulation(Vbg)
+#   bgmodel = get_background_model(Ωbg)
+#   own_model = remove_ghost_cells(bgmodel,gids)
+#   geo1 = DiscreteGeometry(ϕh,own_model)
+#   geo2 = DiscreteGeometry(-ϕh,own_model,name="")
+#   get_geo_params(geo1,geo2,own_model)
+# end
+
+# gids = get_cell_gids(model)
+# contribs = map(local_views(dΩ),local_views(φh),local_views(gids)) do dΩ,φh,gids
+#   _f = u -> j_iem.F(u,dΩ,alt_get_geo_params(u,get_fe_space(φh),gids)[2]...)
+#   return Gridap.Fields.gradient(_f,φh)
+# end;

scripts/Embedded/embedded_measures.jl

@@ -72,7 +72,10 @@ function get_embedded_measures(s::EmbeddedMeasureCache,args...)
   return s.measures
 end
 
-function get_embedded_triangulations(s::EmbeddedMeasureCache)
+function get_embedded_triangulations(s::EmbeddedMeasureCache,φ)
+  trians, measures = get_geo_params(φ,s.space)
+  s.measures = measures
+  s.trians = trians
   return s.trians
 end
 
@@ -106,8 +109,8 @@ function Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh::Vector,K::Int)
   if K == length(uh)
     # Need to test. My hope is that locally this will work - ghosts will be wrong but
     #  will be thrown away later.
-    contribs = map(local_measures,local_views(uh[end]),local_fields) do dΩ,φh,lf
-      _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,F.get_embedded_dΩ(φh,get_fe_space(φh))...)
+    contribs = map(local_measures,local_fields) do dΩ,lf
+      _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,F.get_embedded_dΩ(u,get_fe_space(lf[K]))...)
       return Gridap.Fields.gradient(_f,lf[K])
     end
   else

src/ChainRules.jl

@@ -139,7 +139,8 @@ function StateParamIntegrandWithMeasure(
   U::FESpace,V_φ::FESpace,U_reg::FESpace,
   assem_U::Assembler,assem_deriv::Assembler
 )
-  φ₀, u₀ = zero(V_φ), zero(U)
+  # TODO: Find alternative to this later, is a fix for grad of embedded FEs
+  φ₀, u₀ = interpolate(x->-sqrt((x[1]-1/2)^2+(x[2]-1/2)^2)+0.2,V_φ), zero(U)
   ∂j∂u_vecdata = collect_cell_vector(U,∇(F,[u₀,φ₀],1))
   ∂j∂φ_vecdata = collect_cell_vector(U_reg,∇(F,[u₀,φ₀],2))
   ∂j∂u_vec = allocate_vector(assem_U,∂j∂u_vecdata)