From ec3177226016c7c185bbd83de04f864c9f7a1336 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Mon, 24 Jun 2024 12:24:58 +1000 Subject: [PATCH] added some testing --- .../Embedded/Examples/2d_thermal_with_ad.jl | 95 ++++++++++++++++ .../Examples/2d_thermal_with_ad_MPI.jl | 101 ++++++++++++++++++ .../gridap_embedded_grad_testing_dist_MWE.jl | 24 ++++- scripts/Embedded/embedded_measures.jl | 9 +- src/ChainRules.jl | 3 +- 5 files changed, 225 insertions(+), 7 deletions(-) create mode 100644 scripts/Embedded/Examples/2d_thermal_with_ad.jl create mode 100644 scripts/Embedded/Examples/2d_thermal_with_ad_MPI.jl diff --git a/scripts/Embedded/Examples/2d_thermal_with_ad.jl b/scripts/Embedded/Examples/2d_thermal_with_ad.jl new file mode 100644 index 00000000..a2c9194b --- /dev/null +++ b/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() \ No newline at end of file diff --git a/scripts/Embedded/Examples/2d_thermal_with_ad_MPI.jl b/scripts/Embedded/Examples/2d_thermal_with_ad_MPI.jl new file mode 100644 index 00000000..090694e4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/scripts/Embedded/MWEs/gridap_embedded_grad_testing_dist_MWE.jl b/scripts/Embedded/MWEs/gridap_embedded_grad_testing_dist_MWE.jl index 4115145a..c5c0e4fa 100644 --- a/scripts/Embedded/MWEs/gridap_embedded_grad_testing_dist_MWE.jl +++ b/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) \ No newline at end of file +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; \ No newline at end of file diff --git a/scripts/Embedded/embedded_measures.jl b/scripts/Embedded/embedded_measures.jl index f2d0f904..5b55de6d 100644 --- a/scripts/Embedded/embedded_measures.jl +++ b/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 diff --git a/src/ChainRules.jl b/src/ChainRules.jl index bdf26b7b..4661c6b7 100644 --- a/src/ChainRules.jl +++ b/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)