diff --git a/scripts/Embedded/Examples/2d_thermal copy.jl b/scripts/Embedded/Examples/2d_thermal copy.jl new file mode 100644 index 00000000..5938f841 --- /dev/null +++ b/scripts/Embedded/Examples/2d_thermal copy.jl @@ -0,0 +1,142 @@ +using Gridap,GridapTopOpt +using Gridap.Adaptivity, Gridap.Geometry +using GridapEmbedded, GridapEmbedded.LevelSetCutters + +using GridapTopOpt: StateParamIntegrandWithMeasure + +path="./results/UnfittedFEM_thermal_compliance_ALM/" +n = 51 +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)) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +el_Δ = get_el_Δ(_model) +h = maximum(el_Δ) +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Ω) + +## Levet-set function space and derivative regularisation space +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) +U_reg = TrialFESpace(V_reg,0) +V_φ = TestFESpace(model,reffe_scalar) + +## Levet-set function +φh = interpolate(x->-cos(4π*x[1])*cos(4π*x[2])-0.2,V_φ) +geo = DiscreteGeometry(φh,model) +cutgeo = cut(model,geo) +Ωs = EmbeddedCollection(model,φh) do cutgeo + Ωin = Triangulation(cutgeo,PHYSICAL) + Γg = GhostSkeleton(cutgeo) + n_Γg = get_normal_vector(Γg) + (; + :Ωin => Ωin, + :dΩin => Measure(Ωin,2*order), + :Γg => Γg, + :dΓg => Measure(Γg,2*order), + :n_Γg => get_normal_vector(Γg) + ) +end + +## Weak form +const γg = 0.1 +a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫((γg*h)*jump(Ωs.n_Γg⋅∇(v))*jump(Ωs.n_Γg⋅∇(u)))Ωs.dΓg +l(v,φ) = ∫(v)dΓ_N + +## Optimisation functionals +J(u,φ) = a(u,u,φ) +Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ + +## Setup solver and FE operators +state_collection = EmbeddedCollection(model,φh) do cutgeo + Ωact = Triangulation(cutgeo,ACTIVE) + V = TestFESpace(Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,0.0) + + Ωin = Triangulation(cutgeo,PHYSICAL) + Γg = GhostSkeleton(cutgeo) + n_Γg = get_normal_vector(Γg) + dΩin = Measure(Ωin,2*order) + dΓg = Measure(Γg,2*order) + n_Γg = get_normal_vector(Γg) + + dΩact = Measure(Ωact,2) + Ωact_out = Triangulation(cutgeo,ACTIVE_OUT) + dΩact_out = Measure(Ωact_out,2) + + a(u,v,φ) = ∫(∇(v)⋅∇(u))dΩin + ∫((γg*h)*jump(n_Γg⋅∇(v))*jump(n_Γg⋅∇(u)))dΓg + l( v,φ) = ∫(v)dΓ_N + J(u,φ) = a(u,u,φ) + Vol(u,φ) = ∫(1/vol_D)dΩin - ∫(vf/vol_D)dΩact - ∫(vf/vol_D)dΩact_out + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) + (; + :state_map => state_map, + :J => StateParamIntegrandWithMeasure(J,state_map), + :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),(Vol,)) + ) +end +pcfs = EmbeddedPDEConstrainedFunctionals(state_collection) + +evaluate!(pcfs,φh) + +using Gridap.Arrays +uh = u +ttrian = Ω +strian = get_triangulation(uh) + +D = num_cell_dims(strian) +sglue = get_glue(strian,Val(D)) +tglue = get_glue(ttrian,Val(D)) + +scells = Arrays.IdentityVector(Int32(num_cells(strian))) +mcells = extend(scells,sglue.mface_to_tface) +tcells = lazy_map(Reindex(mcells),tglue.tface_to_mface) +collect(tcells) + + +## Evolution Method +evo = CutFEMEvolve(V_φ,Ωs,dΩ,h) +reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=InteriorPenalty(V_φ)) +ls_evo = UnfittedFEEvolution(evo,reinit) +reinit!(ls_evo,φh) + +## 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) + +Ωs(φh) + +evaluate!(pcfs,φh) + +# ## Optimiser +# optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; +# γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +# for (it,uh,φh) in optimiser +# if iszero(it % iter_mod) +# writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +# writevtk(Ωs.Γ,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) +# writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +# writevtk(Ωs.Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)]) \ No newline at end of file diff --git a/scripts/Embedded/Examples/2d_thermal.jl b/scripts/Embedded/Examples/2d_thermal.jl index 44c2efa1..a241effb 100644 --- a/scripts/Embedded/Examples/2d_thermal.jl +++ b/scripts/Embedded/Examples/2d_thermal.jl @@ -1,94 +1,115 @@ -using Pkg; Pkg.activate() - using Gridap,GridapTopOpt -include("../embedded_measures_AD_DISABLED.jl") +using Gridap.Adaptivity, Gridap.Geometry +using GridapEmbedded, GridapEmbedded.LevelSetCutters -function main() - path="./results/UnfittedFEM_thermal_compliance_ALM/" - n = 200 - 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 +using GridapTopOpt: StateParamIntegrandWithMeasure - 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") +path="./results/UnfittedFEM_thermal_compliance_ALM/" +mkpath(path) +n = 51 +order = 1 +γ = 0.1 +γ_reinit = 0.5 +max_steps = floor(Int,order*minimum(n)/10) +vf = 0.4 +α_coeff = 2max_steps*γ +iter_mod = 1 - ## 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Ω) +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +el_Δ = get_el_Δ(_model) +h = maximum(el_Δ) +h_refine = maximum(el_Δ)/2 +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") - ## 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) +## 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Ω) + +## Levet-set function space and derivative regularisation space +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) +U_reg = TrialFESpace(V_reg,0) +V_φ = TestFESpace(model,reffe_scalar) - φh = interpolate(x->-cos(4π*x[1])*cos(4*pi*x[2])/4-0.2/4,V_φ) - embedded_meas = EmbeddedMeasureCache(φh,V_φ) - update_meas(φ) = update_embedded_measures!(φ,embedded_meas) - get_meas(φ) = get_embedded_measures(φ,embedded_meas) +## Levet-set function +φh = interpolate(x->-cos(4π*x[1])*cos(4π*x[2])-0.2,V_φ) +Ωs = EmbeddedCollection(model,φh) do cutgeo + Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL)) + Γg = DifferentiableTriangulation(GhostSkeleton(cutgeo)) + n_Γg = get_normal_vector(Γg) + Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo)) + (; + :Ωin => Ωin, + :dΩin => Measure(Ωin,2*order), + :Γg => Γg, + :dΓg => Measure(Γg,2*order), + :n_Γg => get_normal_vector(Γg), + :Γ => Γ, + :dΓ => Measure(Γ,2*order) + ) +end - ## Weak form - a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-6*∇(u)⋅∇(v))dΩ2 - l(v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(v)dΓ_N +## Weak form +const γg = 0.1 +a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫((γg*h)*jump(Ωs.n_Γg⋅∇(v))*jump(Ωs.n_Γg⋅∇(u)))Ωs.dΓg +l(v,φ) = ∫(v)dΓ_N - ## Optimisation functionals - J(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u))dΩ1 - dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u)*q)dΓ; - Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf/vol_D)dΩ; - dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)dΓ +## Optimisation functionals +J(u,φ) = ∫(∇(u)⋅∇(u))Ωs.dΩin +Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ +dVol(q,u,φ) = ∫(-1/vol_D*q/(norm ∘ (∇(φ))))Ωs.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) - dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),update_meas) - Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),update_meas) - dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(dΩ,dΓ_N),update_meas) +## Setup solver and FE operators +state_collection = EmbeddedCollection(model,φh) do cutgeo + Ωact = Triangulation(cutgeo,ACTIVE) + V = TestFESpace(Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,0.0) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) + (; + :state_map => state_map, + :J => StateParamIntegrandWithMeasure(J,state_map), + :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,]) + ) +end +pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,)) + +## Evolution Method +evo = CutFEMEvolve(V_φ,Ωs,dΩ,h) +reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=InteriorPenalty(V_φ)) +ls_evo = UnfittedFEEvolution(evo,reinit) +# reinit!(ls_evo,φh) - ## Evolution Method - ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) +## Hilbertian extension-regularisation problems +# α = α_coeff*(h_refine/order)^2 +# a_hilb(p,q) =∫(α*∇(p)⋅∇(q) + p*q)dΩ; +α = α_coeff*h_refine +a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; +vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) - ## 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,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem]) - ## 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) +## Forward problem and derivatives +_j,_c,_dJ,_dC = evaluate!(pcfs,φh) +writevtk(Ω,path*"check_dJ",cellfields=["dJ"=>FEFunction(dJ,U_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) - 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) +## Optimiser +optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh;debug=true, + γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +for (it,uh,φh,state) in optimiser + if iszero(it % iter_mod) + writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"velh"=>FEFunction(V_φ,state.vel)]) + writevtk(Ωs.Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(Ωs.Γ)]) end - it = get_history(optimiser).niter; uh = get_state(pcfs) - _Ω1,_,_Γ = get_embedded_triangulations(embedded_meas) - writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) - writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)]) + write_history(path*"/history.txt",optimiser.history) end - -main() \ No newline at end of file +it = get_history(optimiser).niter; uh = get_state(pcfs) +writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +writevtk(Ωs.Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(Ωs.Γ)]) \ No newline at end of file diff --git a/scripts/Embedded/Examples/2d_elastic_compliance.jl b/scripts/Embedded/Examples/deprecated/2d_elastic_compliance.jl similarity index 100% rename from scripts/Embedded/Examples/2d_elastic_compliance.jl rename to scripts/Embedded/Examples/deprecated/2d_elastic_compliance.jl diff --git a/scripts/Embedded/Examples/2d_inverse_homenisation.jl b/scripts/Embedded/Examples/deprecated/2d_inverse_homenisation.jl similarity index 100% rename from scripts/Embedded/Examples/2d_inverse_homenisation.jl rename to scripts/Embedded/Examples/deprecated/2d_inverse_homenisation.jl diff --git a/scripts/Embedded/Examples/deprecated/2d_thermal.jl b/scripts/Embedded/Examples/deprecated/2d_thermal.jl new file mode 100644 index 00000000..64992920 --- /dev/null +++ b/scripts/Embedded/Examples/deprecated/2d_thermal.jl @@ -0,0 +1,94 @@ +using Pkg; Pkg.activate() + +using Gridap,GridapTopOpt +include("../embedded_measures_AD_DISABLED.jl") + +function main() + path="./results/UnfittedFEM_thermal_compliance_ALM/" + n = 200 + 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(φ) = update_embedded_measures!(φ,embedded_meas) + get_meas(φ) = get_embedded_measures(φ,embedded_meas) + + ## Weak form + a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-6*∇(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 + dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u)*q)dΓ; + Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf/vol_D)dΩ; + dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)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) + dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),update_meas) + Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),update_meas) + dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(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,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem]) + + ## 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) + 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) + writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) + writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)]) +end + +main() \ No newline at end of file diff --git a/scripts/Embedded/Examples/2d_thermal_optimised.jl b/scripts/Embedded/Examples/deprecated/2d_thermal_optimised.jl similarity index 97% rename from scripts/Embedded/Examples/2d_thermal_optimised.jl rename to scripts/Embedded/Examples/deprecated/2d_thermal_optimised.jl index 4c5f5f3d..9763f94f 100644 --- a/scripts/Embedded/Examples/2d_thermal_optimised.jl +++ b/scripts/Embedded/Examples/deprecated/2d_thermal_optimised.jl @@ -1,95 +1,95 @@ -using Pkg; Pkg.activate() - -using Gridap,GridapTopOpt,GridapEmbedded -include("../embedded_measures_AD_DISABLED.jl") - -function main() - path="./results/UnfittedFEM_thermal_compliance_ALM_Optimised_test/" - n = 200 - order = 1 - γ = 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 - - 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(initial_lsf(8,0.2),V_φ) - embedded_meas = EmbeddedMeasureCache(φh,V_φ) - update_meas(φ) = update_embedded_measures!(φ,embedded_meas) - get_meas(φ) = get_embedded_measures(φ,embedded_meas) - - ## Weak form - a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-6*∇(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 - dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u)*q)dΓ; - Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf/vol_D)dΩ; - dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)dΓ - - ## IntegrandWithEmbeddedMeasure - a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas) - l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),get_meas) - J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),get_meas) - dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),get_meas) - Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),get_meas) - dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(dΩ,dΓ_N),get_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,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem]) - - ## 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) - 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) - writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) - writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)]) -end - +using Pkg; Pkg.activate() + +using Gridap,GridapTopOpt,GridapEmbedded +include("../embedded_measures_AD_DISABLED.jl") + +function main() + path="./results/UnfittedFEM_thermal_compliance_ALM_Optimised_test/" + n = 200 + order = 1 + γ = 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 + + 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(initial_lsf(8,0.2),V_φ) + embedded_meas = EmbeddedMeasureCache(φh,V_φ) + update_meas(φ) = update_embedded_measures!(φ,embedded_meas) + get_meas(φ) = get_embedded_measures(φ,embedded_meas) + + ## Weak form + a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-6*∇(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 + dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u)*q)dΓ; + Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf/vol_D)dΩ; + dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)dΓ + + ## IntegrandWithEmbeddedMeasure + a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas) + l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),get_meas) + J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),get_meas) + dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),get_meas) + Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),get_meas) + dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(dΩ,dΓ_N),get_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,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem]) + + ## 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) + 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) + writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) + writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)]) +end + main() \ No newline at end of file diff --git a/scripts/Embedded/Examples/2d_thermal_optimised_MPI.jl b/scripts/Embedded/Examples/deprecated/2d_thermal_optimised_MPI.jl similarity index 97% rename from scripts/Embedded/Examples/2d_thermal_optimised_MPI.jl rename to scripts/Embedded/Examples/deprecated/2d_thermal_optimised_MPI.jl index 71854bee..925cd163 100644 --- a/scripts/Embedded/Examples/2d_thermal_optimised_MPI.jl +++ b/scripts/Embedded/Examples/deprecated/2d_thermal_optimised_MPI.jl @@ -1,100 +1,100 @@ -using Gridap,GridapTopOpt,GridapEmbedded -using GridapDistributed, GridapPETSc, PartitionedArrays - -include("../embedded_measures_AD_DISABLED.jl") - -function main(parts,distribute) - ranks = distribute(LinearIndices((prod(parts),))) - - path="./results/UnfittedFEM_thermal_compliance_ALM_Optimised_MPI_test/" - n = 200 - order = 1 - γ = 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(8,0.2),V_φ) - embedded_meas = EmbeddedMeasureCache(φh,V_φ) - update_meas(φ) = update_embedded_measures!(φ,embedded_meas) - get_meas(φ) = get_embedded_measures(φ,embedded_meas) - - ## Weak form - a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-6*∇(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 - dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u)*q)dΓ; - Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf)dΩ; - dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)dΓ - - ## IntegrandWithEmbeddedMeasure - a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas) - l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),get_meas) - J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),get_meas) - dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),get_meas) - Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),get_meas) - dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(dΩ,dΓ_N),get_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,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem]) - - ## 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) - 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) - 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) +using Gridap,GridapTopOpt,GridapEmbedded +using GridapDistributed, GridapPETSc, PartitionedArrays + +include("../embedded_measures_AD_DISABLED.jl") + +function main(parts,distribute) + ranks = distribute(LinearIndices((prod(parts),))) + + path="./results/UnfittedFEM_thermal_compliance_ALM_Optimised_MPI_test/" + n = 200 + order = 1 + γ = 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(8,0.2),V_φ) + embedded_meas = EmbeddedMeasureCache(φh,V_φ) + update_meas(φ) = update_embedded_measures!(φ,embedded_meas) + get_meas(φ) = get_embedded_measures(φ,embedded_meas) + + ## Weak form + a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-6*∇(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 + dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u)*q)dΓ; + Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf)dΩ; + dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)dΓ + + ## IntegrandWithEmbeddedMeasure + a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas) + l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),get_meas) + J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),get_meas) + dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),get_meas) + Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),get_meas) + dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(dΩ,dΓ_N),get_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,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem]) + + ## 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) + 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) + 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 \ No newline at end of file diff --git a/scripts/Embedded/Examples/2d_thermal_with_ad.jl b/scripts/Embedded/Examples/deprecated/2d_thermal_with_ad.jl similarity index 97% rename from scripts/Embedded/Examples/2d_thermal_with_ad.jl rename to scripts/Embedded/Examples/deprecated/2d_thermal_with_ad.jl index a2c9194b..77a61706 100644 --- a/scripts/Embedded/Examples/2d_thermal_with_ad.jl +++ b/scripts/Embedded/Examples/deprecated/2d_thermal_with_ad.jl @@ -1,95 +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 - +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() \ No newline at end of file diff --git a/scripts/Embedded/Examples/2d_thermal_with_ad_MPI.jl b/scripts/Embedded/Examples/deprecated/2d_thermal_with_ad_MPI.jl similarity index 97% rename from scripts/Embedded/Examples/2d_thermal_with_ad_MPI.jl rename to scripts/Embedded/Examples/deprecated/2d_thermal_with_ad_MPI.jl index 090694e4..a17cfe99 100644 --- a/scripts/Embedded/Examples/2d_thermal_with_ad_MPI.jl +++ b/scripts/Embedded/Examples/deprecated/2d_thermal_with_ad_MPI.jl @@ -1,101 +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) +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 \ No newline at end of file diff --git a/scripts/Embedded/Examples/test.jl b/scripts/Embedded/Examples/test.jl new file mode 100644 index 00000000..4f608202 --- /dev/null +++ b/scripts/Embedded/Examples/test.jl @@ -0,0 +1,148 @@ +using Gridap,GridapTopOpt +using Gridap.Adaptivity, Gridap.Geometry +using GridapEmbedded, GridapEmbedded.LevelSetCutters + +using GridapTopOpt: StateParamIntegrandWithMeasure + +n = 51 +order = 1 + +_model = CartesianDiscreteModel((0,1,0,1),(n,n)) +base_model = UnstructuredDiscreteModel(_model) +ref_model = refine(base_model, refinement_method = "barycentric") +model = ref_model.model +el_Δ = get_el_Δ(_model) +h = maximum(el_Δ) +h_refine = maximum(el_Δ)/2 +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Ω) + +## Levet-set function space and derivative regularisation space +reffe_scalar = ReferenceFE(lagrangian,Float64,order) +V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) +U_reg = TrialFESpace(V_reg,0) +V_φ = TestFESpace(model,reffe_scalar) + +## Levet-set function +φh = interpolate(x->-cos(4π*x[1])*cos(4π*x[2])-0.2,V_φ) +geo = DiscreteGeometry(φh,model) +cutgeo = cut(model,geo) + +x_φ = get_free_dof_values(φh) +idx = findall(isapprox(0.0;atol=10^-10),x_φ) +@assert isempty(idx) + +Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL)) +Γg = DifferentiableTriangulation(GhostSkeleton(cutgeo)) +n_Γg = get_normal_vector(Γg) + +dΩin = Measure(Ωin,2*order) +dΓg = Measure(Γg,2*order) +n_Γg = get_normal_vector(Γg) + +a(u,v,φ) = ∫(∇(v)⋅∇(u))dΩin #+ ∫((γg*h)*jump(n_Γg⋅∇(v))*jump(n_Γg⋅∇(u)))dΓg +l(v,φ) = ∫(v)dΓ_N + +Ωact = Triangulation(cutgeo,ACTIVE) +V = TestFESpace(Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"]) +U = TrialFESpace(V,0.0) + +uhd = zero(U) +# uhd = interpolate(x->1,U) + +# ∇(a,[uhd,uhd,φh],3) +∇(φ->a(uhd,uhd,φ),φh) + + +function generate_model(D,n) + domain = (D==2) ? (0,1,0,1) : (0,1,0,1,0,1) + cell_partition = (D==2) ? (n,n) : (n,n,n) + base_model = UnstructuredDiscreteModel((CartesianDiscreteModel(domain,cell_partition))) + ref_model = refine(base_model, refinement_method = "barycentric") + model = ref_model.model + return model +end +function level_set(shape::Symbol;N=4) + if shape == :square + x -> max(abs(x[1]-0.5),abs(x[2]-0.5))-0.25 # Square + elseif shape == :corner_2d + x -> ((x[1]-0.5)^N+(x[2]-0.5)^N)^(1/N)-0.25 # Curved corner + elseif shape == :diamond + x -> abs(x[1]-0.5)+abs(x[2]-0.5)-0.25-0/n/10 # Diamond + elseif shape == :circle + x -> sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.5223 # Circle + elseif shape == :square_prism + x -> max(abs(x[1]-0.5),abs(x[2]-0.5),abs(x[3]-0.5))-0.25 # Square prism + elseif shape == :corner_3d + x -> ((x[1]-0.5)^N+(x[2]-0.5)^N+(x[3]-0.5)^N)^(1/N)-0.25 # Curved corner + elseif shape == :diamond_prism + x -> abs(x[1]-0.5)+abs(x[2]-0.5)+abs(x[3]-0.5)-0.25-0/n/10 # Diamond prism + elseif shape == :sphere + x -> sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.53 # Sphere + elseif shape == :regular_2d + x -> cos(2π*x[1])*cos(2π*x[2])-0.11 # "Regular" LSF + elseif shape == :regular_3d + x -> cos(2π*x[1])*cos(2π*x[2])*cos(2π*x[3])-0.11 # "Regular" LSF + end +end + + +D = 2 +n = 10 +model = generate_model(D,n) +φ = level_set(:circle) +f = x -> 1.0 +order = 1 +reffe = ReferenceFE(lagrangian,Float64,order) +V_φ = TestFESpace(model,reffe) + +φh = interpolate(φ,V_φ) +fh = interpolate(f,V_φ) + +# Correction if level set is on top of a node +x_φ = get_free_dof_values(φh) +idx = findall(isapprox(0.0;atol=10^-10),x_φ) +!isempty(idx) && @info "Correcting level values!" +x_φ[idx] .+= 10*eps(eltype(x_φ)) + +geo = DiscreteGeometry(φh,model) +cutgeo = cut(model,geo) + +# A.1) Volume integral + +Ω = Triangulation(cutgeo,PHYSICAL_IN) +Ω_AD = DifferentiableTriangulation(Ω) +dΩ = Measure(Ω_AD,2*order) + +Γ = EmbeddedBoundary(cutgeo) +dΓ = Measure(Γ,2*order) + +J_bulk(φ) = ∫(fh)dΩ +dJ_bulk_AD = gradient(J_bulk,φh) +dJ_bulk_AD_vec = assemble_vector(dJ_bulk_AD,V_φ) + +dJ_bulk_exact(q) = ∫(-fh*q/(norm ∘ (∇(φh))))dΓ +dJ_bulk_exact_vec = assemble_vector(dJ_bulk_exact,V_φ) + +@test norm(dJ_bulk_AD_vec - dJ_bulk_exact_vec) < 1e-10 + +# A.2) Volume integral + +g(fh) = ∇(fh)⋅∇(fh) +J_bulk(φ) = ∫(g(fh))dΩ +dJ_bulk_AD = gradient(J_bulk,φh) +dJ_bulk_AD_vec = assemble_vector(dJ_bulk_AD,V_φ) + +dJ_bulk_exact(q) = ∫(-g(fh)*q/(norm ∘ (∇(φh))))dΓ +dJ_bulk_exact_vec = assemble_vector(dJ_bulk_exact,V_φ) + +@test norm(dJ_bulk_AD_vec - dJ_bulk_exact_vec) < 1e-10 \ No newline at end of file diff --git a/src/ChainRules.jl b/src/ChainRules.jl index bcedbabc..13493f41 100644 --- a/src/ChainRules.jl +++ b/src/ChainRules.jl @@ -84,6 +84,7 @@ function StateParamIntegrandWithMeasure( assem_U::Assembler,assem_deriv::Assembler ) φ₀, u₀ = interpolate(x->-sqrt((x[1]-1/2)^2+(x[2]-1/2)^2)+0.2,V_φ), zero(U) + # TODO: Can we make F a dummy functional? ∂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) @@ -383,7 +384,7 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap Optional arguments enable specification of assemblers and linear solvers. """ - function AffineFEStateMap( + function AffineFEStateMap( biform::Function,liform::Function, U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(U,V), @@ -407,13 +408,13 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap K, b = get_matrix(op), get_vector(op) x = allocate_in_domain(K); fill!(x,zero(eltype(x))) ns = numerical_setup(symbolic_setup(ls,K),K) - fwd_caches = (ns,K,b,x,uhd,assem_U) + fwd_caches = (ns,K,b,x,uhd,assem_U,ls) ## Adjoint cache adjoint_K = assemble_matrix((u,v)->biform(v,u,φh),assem_adjoint,V,U) adjoint_x = allocate_in_domain(adjoint_K); fill!(adjoint_x,zero(eltype(adjoint_x))) adjoint_ns = numerical_setup(symbolic_setup(adjoint_ls,adjoint_K),adjoint_K) - adj_caches = (adjoint_ns,adjoint_K,adjoint_x,assem_adjoint) + adj_caches = (adjoint_ns,adjoint_K,adjoint_x,assem_adjoint,adjoint_ls) A,B,C = typeof(biform), typeof(liform), typeof(spaces) D,E,F = typeof(plb_caches),typeof(fwd_caches), typeof(adj_caches) @@ -429,7 +430,7 @@ get_assemblers(m::AffineFEStateMap) = (m.fwd_caches[6],m.plb_caches[2],m.adj_cac function forward_solve!(φ_to_u::AffineFEStateMap,φh) biform, liform = φ_to_u.biform, φ_to_u.liform U, V, _, _ = φ_to_u.spaces - ns, K, b, x, uhd, assem_U = φ_to_u.fwd_caches + ns, K, b, x, uhd, assem_U, _ = φ_to_u.fwd_caches a_fwd(u,v) = biform(u,v,φh) l_fwd(v) = liform(v,φh) @@ -445,7 +446,7 @@ function dRdφ(φ_to_u::AffineFEStateMap,uh,vh,φh) end function update_adjoint_caches!(φ_to_u::AffineFEStateMap,uh,φh) - adjoint_ns, adjoint_K, _, assem_adjoint = φ_to_u.adj_caches + adjoint_ns, adjoint_K, _, assem_adjoint, _ = φ_to_u.adj_caches U, V, _, _ = φ_to_u.spaces assemble_matrix!((u,v) -> φ_to_u.biform(v,u,φh),adjoint_K,assem_adjoint,V,U) numerical_setup!(adjoint_ns,adjoint_K) @@ -453,7 +454,7 @@ function update_adjoint_caches!(φ_to_u::AffineFEStateMap,uh,φh) end function adjoint_solve!(φ_to_u::AffineFEStateMap,du::AbstractVector) - adjoint_ns, _, adjoint_x, _ = φ_to_u.adj_caches + adjoint_ns, _, adjoint_x, _, _ = φ_to_u.adj_caches solve!(adjoint_x,adjoint_ns,du) return adjoint_x end @@ -526,7 +527,8 @@ struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap adjoint_K = assemble_adjoint_matrix(_jac_adj,assem_adjoint,U,V) adjoint_x = allocate_in_domain(adjoint_K); fill!(adjoint_x,zero(eltype(adjoint_x))) adjoint_ns = numerical_setup(symbolic_setup(adjoint_ls,adjoint_K),adjoint_K) - adj_caches = (adjoint_ns,adjoint_K,adjoint_x,assem_adjoint) + adj_caches = (adjoint_ns,adjoint_K,adjoint_x,assem_adjoint,adjoint_ls) + A, B, C = typeof(res), typeof(jac), typeof(spaces) D, E, F = typeof(plb_caches), typeof(fwd_caches), typeof(adj_caches) return new{A,B,C,D,E,F}(res,jac,spaces,plb_caches,fwd_caches,adj_caches) @@ -561,7 +563,7 @@ function dRdφ(φ_to_u::NonlinearFEStateMap,uh,vh,φh) end function update_adjoint_caches!(φ_to_u::NonlinearFEStateMap,uh,φh) - adjoint_ns, adjoint_K, _, assem_adjoint = φ_to_u.adj_caches + adjoint_ns, adjoint_K, _, assem_adjoint, _ = φ_to_u.adj_caches U, V, _, _ = φ_to_u.spaces jac(du,v) = φ_to_u.jac(uh,du,v,φh) assemble_adjoint_matrix!(jac,adjoint_K,assem_adjoint,U,V) @@ -570,7 +572,7 @@ function update_adjoint_caches!(φ_to_u::NonlinearFEStateMap,uh,φh) end function adjoint_solve!(φ_to_u::NonlinearFEStateMap,du::AbstractVector) - adjoint_ns, _, adjoint_x, _ = φ_to_u.adj_caches + adjoint_ns, _, adjoint_x, _, _ = φ_to_u.adj_caches solve!(adjoint_x,adjoint_ns,du) return adjoint_x end @@ -660,13 +662,13 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap b0 = allocate_in_range(K); fill!(b0,zero(eltype(b0))) x = repeated_allocate_in_domain(nblocks,K); fill!(x,zero(eltype(x))) ns = numerical_setup(symbolic_setup(ls,K),K) - fwd_caches = (ns,K,b,x,uhd,assem_U0,b0,assem_U) + fwd_caches = (ns,K,b,x,uhd,assem_U0,b0,assem_U,ls) ## Adjoint cache adjoint_K = assemble_matrix((u,v)->biform(v,u,φh),assem_adjoint,V0,U0) adjoint_x = repeated_allocate_in_domain(nblocks,adjoint_K); fill!(adjoint_x,zero(eltype(adjoint_x))) adjoint_ns = numerical_setup(symbolic_setup(adjoint_ls,adjoint_K),adjoint_K) - adj_caches = (adjoint_ns,adjoint_K,adjoint_x,assem_adjoint) + adj_caches = (adjoint_ns,adjoint_K,adjoint_x,assem_adjoint,adjoint_ls) A,B,C,D = typeof(biform), typeof(liforms), typeof(spaces), typeof(spaces_0) E,F,G = typeof(plb_caches), typeof(fwd_caches), typeof(adj_caches) @@ -746,7 +748,7 @@ get_assemblers(m::RepeatingAffineFEStateMap) = (m.fwd_caches[8],m.plb_caches[2], function forward_solve!(φ_to_u::RepeatingAffineFEStateMap,φh) biform, liforms = φ_to_u.biform, φ_to_u.liform U0, V0 = φ_to_u.spaces_0 - ns, K, b, x, uhd, assem_U0, b0, _ = φ_to_u.fwd_caches + ns, K, b, x, uhd, assem_U0, b0, _, _ = φ_to_u.fwd_caches a_fwd(u,v) = biform(u,v,φh) assemble_matrix!(a_fwd,K,assem_U0,U0,V0) @@ -780,7 +782,7 @@ function dRdφ(φ_to_u::RepeatingAffineFEStateMap,uh,vh,φh) end function update_adjoint_caches!(φ_to_u::RepeatingAffineFEStateMap,uh,φh) - adjoint_ns, adjoint_K, _, assem_adjoint = φ_to_u.adj_caches + adjoint_ns, adjoint_K, _, assem_adjoint, _ = φ_to_u.adj_caches U0, V0 = φ_to_u.spaces_0 assemble_matrix!((u,v) -> φ_to_u.biform(v,u,φh),adjoint_K,assem_adjoint,V0,U0) numerical_setup!(adjoint_ns,adjoint_K) @@ -788,7 +790,7 @@ function update_adjoint_caches!(φ_to_u::RepeatingAffineFEStateMap,uh,φh) end function adjoint_solve!(φ_to_u::RepeatingAffineFEStateMap,du::AbstractBlockVector) - adjoint_ns, _, adjoint_x, _ = φ_to_u.adj_caches + adjoint_ns, _, adjoint_x, _, _ = φ_to_u.adj_caches U0, V0 = φ_to_u.spaces_0 map(repeated_blocks(U0,adjoint_x),repeated_blocks(U0,du)) do xi, dui solve!(xi,adjoint_ns,dui) @@ -796,8 +798,10 @@ function adjoint_solve!(φ_to_u::RepeatingAffineFEStateMap,du::AbstractBlockVect return adjoint_x end +abstract type AbstractPDEConstrainedFunctionals{N} end + """ - struct PDEConstrainedFunctionals{N,A} + struct PDEConstrainedFunctionals{N,A} <: AbstractPDEConstrainedFunctionals{N} An object that computes the objective, constraints, and their derivatives. @@ -847,7 +851,7 @@ The gradient is then ``\\frac{\\partial F}{\\partial \\varphi} = - If `analytic_dJ = nothing` automatic differentiation will be used. - If `analytic_dC[i] = nothing` automatic differentiation will be used for `C[i]`. """ -struct PDEConstrainedFunctionals{N,A} +struct PDEConstrainedFunctionals{N,A} <: AbstractPDEConstrainedFunctionals{N} J C dJ @@ -895,24 +899,25 @@ Create an instance of `PDEConstrainedFunctionals` when the problem has no constr PDEConstrainedFunctionals(J,state_map::AbstractFEStateMap;analytic_dJ=nothing) = PDEConstrainedFunctionals(J,typeof(J)[],state_map;analytic_dJ = analytic_dJ,analytic_dC = Nothing[]) -get_state(m::PDEConstrainedFunctionals) = get_state(m.state_map) +get_state_map(m::PDEConstrainedFunctionals) = m.state_map +get_state(m::PDEConstrainedFunctionals) = get_state(get_state_map(m)) """ evaluate_functionals!(pcf::PDEConstrainedFunctionals,φh) Evaluate the objective and constraints at `φh`. """ -function evaluate_functionals!(pcf::PDEConstrainedFunctionals,φh) - u = pcf.state_map(φh) - U = get_trial_space(pcf.state_map) +function evaluate_functionals!(pcf::PDEConstrainedFunctionals,φh;kwargs...) + u = get_state_map(pcf)(φh) + U = get_trial_space(get_state_map(pcf)) uh = FEFunction(U,u) return pcf.J(uh,φh), map(Ci->Ci(uh,φh),pcf.C) end -function evaluate_functionals!(pcf::PDEConstrainedFunctionals,φ::AbstractVector) - V_φ = get_aux_space(pcf.state_map) +function evaluate_functionals!(pcf::PDEConstrainedFunctionals,φ::AbstractVector;kwargs...) + V_φ = get_aux_space(get_state_map(pcf)) φh = FEFunction(V_φ,φ) - return evaluate_functionals!(pcf,φh) + return evaluate_functionals!(pcf,φh;kwargs...) end """ @@ -920,15 +925,15 @@ end Evaluate the derivatives of the objective and constraints at `φh`. """ -function evaluate_derivatives!(pcf::PDEConstrainedFunctionals,φh) +function evaluate_derivatives!(pcf::PDEConstrainedFunctionals,φh;kwargs...) _,_,dJ,dC = evaluate!(pcf,φh) return dJ,dC end -function evaluate_derivatives!(pcf::PDEConstrainedFunctionals,φ::AbstractVector) - V_φ = get_aux_space(pcf.state_map) +function evaluate_derivatives!(pcf::PDEConstrainedFunctionals,φ::AbstractVector;kwargs...) + V_φ = get_aux_space(get_state_map(pcf)) φh = FEFunction(V_φ,φ) - return evaluate_derivatives!(pcf,φh) + return evaluate_derivatives!(pcf,φh;kwargs...) end """ @@ -937,17 +942,179 @@ end Evaluate the objective and constraints, and their derivatives at `φh`. """ -function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh) +function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh;kwargs...) J, C, dJ, dC = pcf.J,pcf.C,pcf.dJ,pcf.dC analytic_dJ = pcf.analytic_dJ analytic_dC = pcf.analytic_dC - U = get_trial_space(pcf.state_map) + U = get_trial_space(get_state_map(pcf)) + + U_reg = get_deriv_space(get_state_map(pcf)) + deriv_assem = get_deriv_assembler(get_state_map(pcf)) + + ## Foward problem + u, u_pullback = rrule(get_state_map(pcf),φh) + uh = FEFunction(U,u) + + function ∇!(F::StateParamIntegrandWithMeasure,dF,::Nothing) + # Automatic differentation + j_val, j_pullback = rrule(F,uh,φh) # Compute functional and pull back + _, dFdu, dFdφ = j_pullback(1) # Compute dFdu, dFdφ + _, dφ_adj = u_pullback(dFdu) # Compute -dFdu*dudφ via adjoint + copy!(dF,dφ_adj) + dF .+= dFdφ + return j_val + end + function ∇!(F::StateParamIntegrandWithMeasure,dF,dF_analytic::Function) + # Analytic shape derivative + j_val = F(uh,φh) + _dF(q) = dF_analytic(q,uh,φh) + assemble_vector!(_dF,dF,deriv_assem,U_reg) + return j_val + end + j = ∇!(J,dJ,analytic_dJ) + c = map(∇!,C,dC,analytic_dC) + + return j,c,dJ,dC +end + +function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φ::AbstractVector;kwargs...) + V_φ = get_aux_space(get_state_map(pcf)) + φh = FEFunction(V_φ,φ) + return evaluate!(pcf,φh;kwargs...) +end + +""" + mutable struct EmbeddedPDEConstrainedFunctionals{N} <: AbstractPDEConstrainedFunctionals{N} + +A mutable version of `PDEConstrainedFunctionals` that allows `state_map` to be +updated given new FE spaces for the forward problem. This is currently required +for body-fitted mesh methods and unfitted methods. +""" +struct EmbeddedPDEConstrainedFunctionals{N,T} <: AbstractPDEConstrainedFunctionals{N} + dJ + dC + analytic_dJ + analytic_dC + embedded_collection + + @doc """ + EmbeddedPDEConstrainedFunctionals(objective::Function,constraints::Vector{<:Function}, + embedded_collection :: EmbeddedCollection;analytic_dJ;analytic_dC) + + Create an instance of `EmbeddedPDEConstrainedFunctionals`. + """ + function EmbeddedPDEConstrainedFunctionals( + embedded_collection :: EmbeddedCollection; + analytic_dJ = nothing, + analytic_dC = nothing) + + @assert Set((:state_map,:J,:C)) == keys(embedded_collection.objects) """ + Expected EmbeddedCollection to have objects ':state_map,:J,:C'. Ensure that you + have updated the collection after adding new recipes. + + You have $(keys(embedded_collection.objects)) + + Note: + - We require that this EmbeddedCollection is seperate to the one used for the + UnfittedEvolution. This is because updating the FEStateMap is more expensive than + cutting and there are instances where evolution and reinitialisation happen + at before recomputing the forward solution. As such, we cut an extra time + to avoid allocating the state map more often then required. + - For problems with no constraints `:C` must at least point to an empty list + """ + # Preallocate + dJ = similar(embedded_collection.J.caches[2]) + dC = map(Ci->similar(Ci.caches[2]),embedded_collection.C) + + N = length(embedded_collection.C) + if analytic_dC isa Nothing + analytic_dC = fill(nothing,length(N)) + end + + T = typeof(embedded_collection.state_map) + return new{N,T}(dJ,dC,analytic_dJ,analytic_dC,embedded_collection) + end +end + +get_state_map(m::EmbeddedPDEConstrainedFunctionals) = m.embedded_collection.state_map +get_state(m::EmbeddedPDEConstrainedFunctionals) = get_state(get_state_map(m)) + +""" + evaluate_functionals!(pcf::EmbeddedPDEConstrainedFunctionals,φh) + +Evaluate the objective and constraints at `φh`. + +!!! warning + Taking `update_space = false` will NOT update the underlying finite element + spaces and assemblers that depend on `φh`. This should only be used + when you are certain that `φh` has not been updated. +""" +function evaluate_functionals!(pcf::EmbeddedPDEConstrainedFunctionals,φh;update_space::Bool=true) + update_space && update_collection!(pcf.embedded_collection,φh) + u = get_state_map(pcf)(φh) + U = get_trial_space(get_state_map(pcf)) + uh = FEFunction(U,u) + J = pcf.embedded_collection.J + C = pcf.embedded_collection.C + return J(uh,φh), map(Ci->Ci(uh,φh),C) +end + +function evaluate_functionals!(pcf::EmbeddedPDEConstrainedFunctionals,φ::AbstractVector;kwargs...) + V_φ = get_aux_space(get_state_map(pcf)) + φh = FEFunction(V_φ,φ) + return evaluate_functionals!(pcf,φh;kwargs...) +end + +""" + evaluate_derivatives!(pcf::EmbeddedPDEConstrainedFunctionals,φh) + +Evaluate the derivatives of the objective and constraints at `φh`. + +!!! warning + Taking `update_space = false` will NOT update the underlying finite element + spaces and assemblers that depend on `φh`. This should only be used + when you are certain that `φh` has not been updated. +""" +function evaluate_derivatives!(pcf::EmbeddedPDEConstrainedFunctionals,φh;update_space::Bool=true) + update_space && update_collection!(pcf.embedded_collection,φh) + _,_,dJ,dC = evaluate!(pcf,φh) + return dJ,dC +end + +function evaluate_derivatives!(pcf::EmbeddedPDEConstrainedFunctionals,φ::AbstractVector;kwargs...) + V_φ = get_aux_space(get_state_map(pcf)) + φh = FEFunction(V_φ,φ) + return evaluate_derivatives!(pcf,φh;kwargs...) +end + +""" + Fields.evaluate!(pcf::EmbeddedPDEConstrainedFunctionals,φh) + +Evaluate the objective and constraints, and their derivatives at +`φh`. + +!!! warning + Taking `update_space = false` will NOT update the underlying finite element + spaces and assemblers that depend on `φh`. This should only be used + when you are certain that `φh` has not been updated. +""" +function Fields.evaluate!(pcf::EmbeddedPDEConstrainedFunctionals,φh;update_space::Bool=true) + update_space && update_collection!(pcf.embedded_collection,φh) + + J = pcf.embedded_collection.J + C = pcf.embedded_collection.C + dJ = pcf.dJ + dC = pcf.dC + analytic_dJ = pcf.analytic_dJ + analytic_dC = pcf.analytic_dC + state_map = get_state_map(pcf) + U = get_trial_space(state_map) - U_reg = get_deriv_space(pcf.state_map) - deriv_assem = get_deriv_assembler(pcf.state_map) + U_reg = get_deriv_space(state_map) + deriv_assem = get_deriv_assembler(state_map) ## Foward problem - u, u_pullback = rrule(pcf.state_map,φh) + u, u_pullback = rrule(state_map,φh) uh = FEFunction(U,u) function ∇!(F::StateParamIntegrandWithMeasure,dF,::Nothing) @@ -972,10 +1139,10 @@ function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh) return j,c,dJ,dC end -function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φ::AbstractVector) - V_φ = get_aux_space(pcf.state_map) +function Fields.evaluate!(pcf::EmbeddedPDEConstrainedFunctionals,φ::AbstractVector;kwargs...) + V_φ = get_aux_space(get_state_map(pcf)) φh = FEFunction(V_φ,φ) - return evaluate!(pcf,φh) + return evaluate!(pcf,φh;kwargs...) end # IO @@ -988,7 +1155,7 @@ function Base.show(io::IO,object::AbstractFEStateMap) print(io,"$(nameof(typeof(object)))") end -function Base.show(io::IO,::MIME"text/plain",f::PDEConstrainedFunctionals) +function Base.show(io::IO,::MIME"text/plain",f::AbstractPDEConstrainedFunctionals{N}) where N print(io,"$(nameof(typeof(f))): - num_constraints: $(length(f.C))") + num_constraints: $N") end diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl index ba410c04..306e0763 100644 --- a/src/GridapTopOpt.jl +++ b/src/GridapTopOpt.jl @@ -44,8 +44,14 @@ import Gridap: solve! include("GridapExtensions.jl") +include("Embedded/Embedded.jl") +export DifferentiableTriangulation +export SubFacetBoundaryTriangulation, SubFacetSkeletonTriangulation +export EmbeddedCollection, update_collection!, add_recipe! + include("ChainRules.jl") export PDEConstrainedFunctionals +export EmbeddedPDEConstrainedFunctionals export AffineFEStateMap export NonlinearFEStateMap export RepeatingAffineFEStateMap @@ -53,11 +59,6 @@ export get_state export evaluate_functionals! export evaluate_derivatives! -include("Embedded/Embedded.jl") -export DifferentiableTriangulation -export SubFacetBoundaryTriangulation, SubFacetSkeletonTriangulation -export EmbeddedCollection, update_collection!, add_recipe! - include("Utilities.jl") export SmoothErsatzMaterialInterpolation export update_labels! diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl index d8710f6d..820a2ca7 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl @@ -121,7 +121,7 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing) _, φhF = data copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) s.cache = state_new - + update_collection!(s.Ωs,φh) # TODO: remove? return φh end @@ -160,7 +160,7 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache) _, φhF = data copy!(get_free_dof_values(φh),get_free_dof_values(φhF)) s.cache = state_updated - + update_collection!(s.Ωs,φh) # TODO: remove? return φh end diff --git a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl index 4bda0a3e..3142f8e4 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl @@ -76,6 +76,7 @@ function solve!(s::StabilisedReinit,φh,nls_cache::Nothing) s.cache = nls_cache # Solve solve!(get_free_dof_values(φh),nls,op,nls_cache) + update_collection!(s.Ωs,φh) # TODO: remove? return φh end @@ -87,6 +88,7 @@ function solve!(s::StabilisedReinit,φh,nls_cache) op = get_algebraic_operator(FEOperator(res,jac,V_φ,V_φ,assembler)) # Solve solve!(get_free_dof_values(φh),nls,op,nls_cache) + update_collection!(s.Ωs,φh) # TODO: remove? return φh end @@ -131,7 +133,6 @@ function get_residual_and_jacobian(s::StabilisedReinit{InteriorPenalty},φh) dΓ = Ωs.dΓ W(u,∇u) = sign(u) * ∇u / (ϵ + norm(∇u)) - V(w) = ca*h*(sqrt ∘ ( w ⋅ w )) a(w,u,v) = ∫(v*(W ∘ (w,∇(w))) ⋅ ∇(u))dΩ_bg + ∫(h^2*jump(∇(u)) ⋅ jump(∇(v)))dΛ + ∫((γd/h)*v*u)dΓ b(w,v) = ∫((sign ∘ w)*v)dΩ_bg res(u,v) = a(u,u,v) - b(u,v) diff --git a/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl b/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl index 47d279a8..7ebbd6d0 100644 --- a/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl +++ b/src/LevelSetEvolution/UnfittedEvolution/UnfittedEvolution.jl @@ -51,9 +51,15 @@ function evolve!(s::UnfittedFEEvolution,φ::T,vel::M,γ,args...) where end function evolve!(s::UnfittedFEEvolution,φ::T,velh,γ,args...) where T<:AbstractVector -V_φ = get_space(s.evolver) -φh = FEFunction(V_φ,φ) -evolve!(s,φh,velh,γ,args...) + V_φ = get_space(s.evolver) + φh = FEFunction(V_φ,φ) + evolve!(s,φh,velh,γ,args...) +end + +function evolve!(s::UnfittedFEEvolution,φh,vel::T,γ,args...) where T<:AbstractVector + V_φ = get_space(s.evolver) + velh = FEFunction(V_φ,vel) + evolve!(s,φh,velh,γ,args...) end function evolve!(s::UnfittedFEEvolution,φh,velh,γ,args...) diff --git a/src/Optimisers/AugmentedLagrangian.jl b/src/Optimisers/AugmentedLagrangian.jl index bf33e91d..d8473b45 100644 --- a/src/Optimisers/AugmentedLagrangian.jl +++ b/src/Optimisers/AugmentedLagrangian.jl @@ -4,11 +4,11 @@ An augmented Lagrangian method based on Nocedal and Wright, 2006 ([link](https://doi.org/10.1007/978-0-387-40065-5)). Note that this method will function as a Lagrangian method if no constraints -are defined in `problem::PDEConstrainedFunctionals`. +are defined in `problem::AbstractPDEConstrainedFunctionals`. # Parameters -- `problem::PDEConstrainedFunctionals`: The objective and constraint setup. +- `problem::AbstractPDEConstrainedFunctionals`: The objective and constraint setup. - `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. @@ -22,7 +22,7 @@ iteration history. By default this uses a mean zero crossing algorithm as implem in ChaosTools. Oscillations checking can be disabled by taking `has_oscillations = (args...) -> false`. """ struct AugmentedLagrangian <: Optimiser - problem :: PDEConstrainedFunctionals + problem :: AbstractPDEConstrainedFunctionals ls_evolver :: LevelSetEvolution vel_ext :: VelocityExtension history :: OptimiserHistory{Float64} @@ -33,7 +33,7 @@ struct AugmentedLagrangian <: Optimiser @doc """ AugmentedLagrangian( - problem :: PDEConstrainedFunctionals{N}, + problem :: AbstractPDEConstrainedFunctionals{N}, ls_evolver :: LevelSetEvolution, vel_ext :: VelocityExtension, φ0; @@ -47,7 +47,7 @@ struct AugmentedLagrangian <: Optimiser # Required - - `problem::PDEConstrainedFunctionals`: The objective and constraint setup. + - `problem::AbstractPDEConstrainedFunctionals`: The objective and constraint setup. - `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. @@ -73,7 +73,7 @@ struct AugmentedLagrangian <: Optimiser - `debug = false`: Debug flag. """ function AugmentedLagrangian( - problem :: PDEConstrainedFunctionals{N}, + problem :: AbstractPDEConstrainedFunctionals{N}, ls_evolver :: LevelSetEvolution, vel_ext :: VelocityExtension, φ0; @@ -191,8 +191,8 @@ function Base.iterate(m::AugmentedLagrangian,state) print_msg(m.history," Oscillations detected, reducing γ to $(γ)\n",color=:yellow) end - U_reg = get_deriv_space(m.problem.state_map) - V_φ = get_aux_space(m.problem.state_map) + U_reg = get_deriv_space(get_state_map(m.problem)) + V_φ = get_aux_space(get_state_map(m.problem)) interpolate!(FEFunction(U_reg,dL),vel,V_φ) evolve!(m.ls_evolver,φh,vel,γ) iszero(it % reinit_mod) && reinit!(m.ls_evolver,φh,γ_reinit) diff --git a/src/Optimisers/HilbertianProjection.jl b/src/Optimisers/HilbertianProjection.jl index 9b576127..7e69513d 100644 --- a/src/Optimisers/HilbertianProjection.jl +++ b/src/Optimisers/HilbertianProjection.jl @@ -138,7 +138,7 @@ A Hilbertian projection method as described by Wegert et al., 2023 # Parameters -- `problem::PDEConstrainedFunctionals{N}`: The objective and constraint setup. +- `problem::AbstractPDEConstrainedFunctionals{N}`: The objective and constraint setup. - `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. @@ -149,7 +149,7 @@ A Hilbertian projection method as described by Wegert et al., 2023 - `params::NamedTuple`: Optimisation parameters. """ struct HilbertianProjection{A} <: Optimiser - problem :: PDEConstrainedFunctionals + problem :: AbstractPDEConstrainedFunctionals ls_evolver :: LevelSetEvolution vel_ext :: VelocityExtension projector :: HilbertianProjectionMap @@ -161,7 +161,7 @@ struct HilbertianProjection{A} <: Optimiser @doc """ HilbertianProjection( - problem :: PDEConstrainedFunctionals{N}, + problem :: AbstractPDEConstrainedFunctionals{N}, ls_evolver :: LevelSetEvolution, vel_ext :: VelocityExtension, φ0; @@ -179,7 +179,7 @@ struct HilbertianProjection{A} <: Optimiser # Required - - `problem::PDEConstrainedFunctionals{N}`: The objective and constraint setup. + - `problem::AbstractPDEConstrainedFunctionals{N}`: The objective and constraint setup. - `ls_evolver::LevelSetEvolution`: Solver for the evolution and reinitisation equations. - `vel_ext::VelocityExtension`: The velocity-extension method for extending shape sensitivities onto the computational domain. @@ -229,7 +229,7 @@ struct HilbertianProjection{A} <: Optimiser This can be set to always be enfored by taking `ls_ξ = 0.0025` and `ls_ξ_reduce_coef = 0.1`. """ function HilbertianProjection( - problem :: PDEConstrainedFunctionals{N}, + problem :: AbstractPDEConstrainedFunctionals{N}, ls_evolver :: LevelSetEvolution, vel_ext :: VelocityExtension, φ0; @@ -349,8 +349,8 @@ function Base.iterate(m::HilbertianProjection,state) end ## Line search - U_reg = get_deriv_space(m.problem.state_map) - V_φ = get_aux_space(m.problem.state_map) + U_reg = get_deriv_space(get_state_map(m.problem)) + V_φ = get_aux_space(get_state_map(m.problem)) interpolate!(FEFunction(U_reg,θ),vel,V_φ) J, C, dJ, dC = _linesearch!(m,state,γ) diff --git a/src/Utilities.jl b/src/Utilities.jl index a3bc5786..1528eb81 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -69,13 +69,13 @@ end """ update_labels!(e::Int,model,f_Γ::Function,name::String) -Given a tag number `e`, a `CartesianDiscreteModel` or `DistributedDiscreteModel` model, -an indicator function `f_Γ`, and a string `name`, label the corresponding vertices, edges, and faces +Given a tag number `e`, a `DiscreteModel` model, an indicator function `f_Γ`, +and a string `name`, label the corresponding vertices, edges, and faces as `name`. Note: `f_Γ` must recieve a Vector and return a Boolean depending on whether it indicates Γ """ -function update_labels!(e::Integer,model::CartesianDiscreteModel,f_Γ::Function,name::String) +function update_labels!(e::Integer,model::DiscreteModel,f_Γ::Function,name::String) mask = mark_nodes(f_Γ,model) _update_labels_locally!(e,model,mask,name) nothing @@ -92,7 +92,7 @@ function update_labels!(e::Integer,model::DistributedDiscreteModel,f_Γ::Functio nothing end -function _update_labels_locally!(e,model::CartesianDiscreteModel{2},mask,name) +function _update_labels_locally!(e,model::DiscreteModel{2},mask,name) topo = get_grid_topology(model) labels = get_face_labeling(model) cell_to_entity = labels.d_to_dface_to_entity[end] @@ -110,7 +110,7 @@ function _update_labels_locally!(e,model::CartesianDiscreteModel{2},mask,name) return cell_to_entity end -function _update_labels_locally!(e,model::CartesianDiscreteModel{3},mask,name) +function _update_labels_locally!(e,model::DiscreteModel{3},mask,name) topo = get_grid_topology(model) labels = get_face_labeling(model) cell_to_entity = labels.d_to_dface_to_entity[end] @@ -164,8 +164,7 @@ initial_lsf(ξ,a;b=0) = x::VectorValue -> -1/4*prod(cos.(get_array(@.(ξ*pi*(x-b """ get_el_Δ(model) -Given a CartesianDiscreteModel or DistributedDiscreteModel that is -uniform, return the element size as a tuple. +Given a CartesianDiscreteModel return the element size as a tuple. """ function get_el_Δ(model::CartesianDiscreteModel) desc = get_cartesian_descriptor(model) diff --git a/test/seq/EmbeddedDifferentiationTests.jl b/test/seq/EmbeddedDifferentiationTests.jl index 2c4dc7a0..a428a722 100644 --- a/test/seq/EmbeddedDifferentiationTests.jl +++ b/test/seq/EmbeddedDifferentiationTests.jl @@ -62,7 +62,7 @@ function main( geo = DiscreteGeometry(φh,model) cutgeo = cut(model,geo) - # A) Volume integral + # A.1) Volume integral Ω = Triangulation(cutgeo,PHYSICAL_IN) Ω_AD = DifferentiableTriangulation(Ω) @@ -80,7 +80,19 @@ function main( @test norm(dJ_bulk_AD_vec - dJ_bulk_exact_vec) < 1e-10 - # B) Facet integral + # A.2) Volume integral + + g(fh) = ∇(fh)⋅∇(fh) + J_bulk2(φ) = ∫(g(fh))dΩ + dJ_bulk_AD2 = gradient(J_bulk2,φh) + dJ_bulk_AD_vec2 = assemble_vector(dJ_bulk_AD2,V_φ) + + dJ_bulk_exact2(q) = ∫(-g(fh)*q/(norm ∘ (∇(φh))))dΓ + dJ_bulk_exact_vec2 = assemble_vector(dJ_bulk_exact2,V_φ) + + @test norm(dJ_bulk_AD_vec2 - dJ_bulk_exact_vec2) < 1e-10 + + # B.1) Facet integral Γ = EmbeddedBoundary(cutgeo) Γ_AD = DifferentiableTriangulation(Γ) @@ -113,6 +125,19 @@ function main( @test norm(dJ_int_AD_vec - dJ_int_exact_vec) < 1e-10 + # B.2) Facet integral + + J_int2(φ) = ∫(g(fh))dΓ_AD + dJ_int_AD2 = gradient(J_int2,φh) + dJ_int_AD_vec2 = assemble_vector(dJ_int_AD2,V_φ) + + dJ_int_exact2(w) = ∫((-n_Γ⋅∇(g(fh)))*w/(norm ∘ (∇(φh))))dΓ + + ∫(-n_S_Λ ⋅ (jump(g(fh)*m_k_Λ) * mean(w) / ∇ˢφ_Λ))dΛ + + ∫(-n_S_Σ ⋅ (g(fh)*m_k_Σ * w / ∇ˢφ_Σ))dΣ + dJ_int_exact_vec2 = assemble_vector(dJ_int_exact2,V_φ) + + @test norm(dJ_int_AD_vec2 - dJ_int_exact_vec2) < 1e-10 + if vtk path = "results/$(name)" Ω_bg = Triangulation(model) @@ -179,4 +204,11 @@ model = generate_model(D,n) f = x -> 1.0 main(model,φ,f;vtk=true) +D = 2 +n = 10 +model = generate_model(D,n) +φ = level_set(:circle) +f = x -> x[1]+x[2] +main(model,φ,f;vtk=true) + end \ No newline at end of file