From 28ff343201fea5bef07c3e3031400f6d10c175e5 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 12 Mar 2024 12:57:20 +1100 Subject: [PATCH 01/14] Added mem leak script --- scripts/_dev/mem_leak.jl | 114 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 scripts/_dev/mem_leak.jl diff --git a/scripts/_dev/mem_leak.jl b/scripts/_dev/mem_leak.jl new file mode 100644 index 00000000..f4fc3270 --- /dev/null +++ b/scripts/_dev/mem_leak.jl @@ -0,0 +1,114 @@ +using Gridap, LevelSetTopOpt + +# FE parameters +order = 1 # Finite element order +xmax=ymax=zmax=1.0 # Domain size +dom = (0,xmax,0,ymax,0,zmax) # Bounding domain +el_size = (40,40,40) # Mesh partition size +prop_Γ_N = 0.4 # Γ_N size parameter +prop_Γ_D = 0.2 # Γ_D size parameter +f_Γ_N(x) = (x[1] ≈ xmax) && # Γ_N indicator function + (ymax/2-ymax*prop_Γ_N/4 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/4 + eps()) && + (zmax/2-zmax*prop_Γ_N/4 - eps() <= x[3] <= zmax/2+zmax*prop_Γ_N/4 + eps()) +f_Γ_D(x) = (x[1] ≈ 0.0) && # Γ_D indicator function + (x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps()) && + (x[3] <= zmax*prop_Γ_D + eps() || x[3] >= zmax-zmax*prop_Γ_D - eps()) +# FD parameters +γ = 0.1 # HJ equation time step coefficient +γ_reinit = 0.5 # Reinit. equation time step coefficient +max_steps = floor(Int,minimum(el_size)/3) # Max steps for advection +tol = 1/(2order^2)*prod(inv,minimum(el_size)) # Advection tolerance +# Problem parameters +κ = 1 # Diffusivity +g = 1 # Heat flow in +vf = 0.4 # Volume fraction constraint +lsf_func = initial_lsf(4,0.2) # Initial level set function +iter_mod = 10 # Output VTK files every 10th iteration + +# Model +model = CartesianDiscreteModel(dom,el_size); +update_labels!(1,model,f_Γ_D,"Gamma_D") +update_labels!(2,model,f_Γ_N,"Gamma_N") +# Triangulation and measures +Ω = Triangulation(model) +Γ_N = BoundaryTriangulation(model,tags="Gamma_N") +dΩ = Measure(Ω,2*order) +dΓ_N = Measure(Γ_N,2*order) +# Spaces +reffe = ReferenceFE(lagrangian,Float64,order) +V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) +U = TrialFESpace(V,0.0) +V_φ = TestFESpace(model,reffe) +V_reg = TestFESpace(model,reffe;dirichlet_tags=["Gamma_N"]) +U_reg = TrialFESpace(V_reg,0) +# Level set and interpolator +φh = interpolate(lsf_func,V_φ) +interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model))) +I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ +# Weak formulation +a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ +l(v,φ,dΩ,dΓ_N) = ∫(g*v)dΓ_N +state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) +# Objective and constraints +J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ +dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ +vol_D = sum(∫(1)dΩ) +C(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ +dC(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ +pcfs = PDEConstrainedFunctionals(J,[C],state_map,analytic_dJ=dJ,analytic_dC=[dC]) +# Velocity extension +α = 4*maximum(get_el_Δ(model)) +a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ +vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) +# Finite difference scheme +scheme = FirstOrderStencil(length(el_size),Float64) +stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) +# Optimiser +optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true,constraint_names=[:Vol]) +# Solve +for i = 1:10 + LevelSetTopOpt.rrule(state_map,φh) +end + +evaluate!(pcfs,φh); +uh = get_state(state_map) + +GC.gc() +Sys.free_memory() / 2^20 + +for i = 1:10 + begin + u = LevelSetTopOpt.forward_solve!(state_map,φh); + uh = FEFunction(LevelSetTopOpt.get_trial_space(state_map),u); + # LevelSetTopOpt.update_adjoint_caches!(state_map,uh,φh); + end; +end; + +Sys.free_memory() / 2^20 + +GC.gc() +Sys.free_memory() / 2^20 + +for i = 1:10 + begin + # u = LevelSetTopOpt.forward_solve!(state_map,φh); + # uh = FEFunction(LevelSetTopOpt.get_trial_space(state_map),u); + + LevelSetTopOpt.update_adjoint_caches!(state_map,uh,φh); + end; +end; + +Sys.free_memory() / 2^20 + +GC.gc() +Sys.free_memory() / 2^20 + +for i = 1:10 + begin + u = LevelSetTopOpt.forward_solve!(state_map,φh); + uh = FEFunction(LevelSetTopOpt.get_trial_space(state_map),u); + LevelSetTopOpt.update_adjoint_caches!(state_map,uh,φh); + end; +end; + +Sys.free_memory() / 2^20 From ab8127a05a5989eda2d29a4fa5a10b105b7cfdf8 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 12 Mar 2024 12:58:22 +1100 Subject: [PATCH 02/14] Minor --- scripts/_dev/mem_leak.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/_dev/mem_leak.jl b/scripts/_dev/mem_leak.jl index f4fc3270..af9a2b28 100644 --- a/scripts/_dev/mem_leak.jl +++ b/scripts/_dev/mem_leak.jl @@ -4,7 +4,7 @@ using Gridap, LevelSetTopOpt order = 1 # Finite element order xmax=ymax=zmax=1.0 # Domain size dom = (0,xmax,0,ymax,0,zmax) # Bounding domain -el_size = (40,40,40) # Mesh partition size +el_size = (10,10,10) # Mesh partition size prop_Γ_N = 0.4 # Γ_N size parameter prop_Γ_D = 0.2 # Γ_D size parameter f_Γ_N(x) = (x[1] ≈ xmax) && # Γ_N indicator function @@ -17,7 +17,7 @@ f_Γ_D(x) = (x[1] ≈ 0.0) && # Γ_D indicator func γ = 0.1 # HJ equation time step coefficient γ_reinit = 0.5 # Reinit. equation time step coefficient max_steps = floor(Int,minimum(el_size)/3) # Max steps for advection -tol = 1/(2order^2)*prod(inv,minimum(el_size)) # Advection tolerance +tol = 1/(2*order^2)*prod(inv,minimum(el_size)) # Advection tolerance # Problem parameters κ = 1 # Diffusivity g = 1 # Heat flow in From 7907edb9fb117f35d5f243cf62fe1368e72c056c Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 12 Mar 2024 13:03:09 +1100 Subject: [PATCH 03/14] Better gitignores --- .gitignore | 2 -- results/.gitignore | 2 ++ results/empty | 1 - 3 files changed, 2 insertions(+), 3 deletions(-) create mode 100644 results/.gitignore delete mode 100644 results/empty diff --git a/.gitignore b/.gitignore index b8aba4f7..d43b6663 100644 --- a/.gitignore +++ b/.gitignore @@ -19,8 +19,6 @@ deps/src/ # Build artifacts for creating documentation generated by the Documenter package docs/build/ docs/site/ -Results/* -results/* # File generated by Pkg, the package manager, based on a corresponding Project.toml # It records a fixed state of all packages used by the project. As such, it should not be diff --git a/results/.gitignore b/results/.gitignore new file mode 100644 index 00000000..c96a04f0 --- /dev/null +++ b/results/.gitignore @@ -0,0 +1,2 @@ +* +!.gitignore \ No newline at end of file diff --git a/results/empty b/results/empty deleted file mode 100644 index 0519ecba..00000000 --- a/results/empty +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file From fa554e9a6c893f88456ff9de3e8294af8fe0aa16 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Tue, 12 Mar 2024 13:54:38 +1000 Subject: [PATCH 04/14] Typo fix --- scripts/MPI/inverse_homenisation_ALM.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/MPI/inverse_homenisation_ALM.jl b/scripts/MPI/inverse_homenisation_ALM.jl index f08691d0..ed9c279f 100644 --- a/scripts/MPI/inverse_homenisation_ALM.jl +++ b/scripts/MPI/inverse_homenisation_ALM.jl @@ -44,7 +44,7 @@ function main(mesh_partition,distribute,el_size) reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) reffe_scalar = ReferenceFE(lagrangian,Float64,order) V = TestFESpace(model,reffe;dirichlet_tags=["origin"]) - U = TrialFESpace(_V,VectorValue(0.0,0.0)) + U = TrialFESpace(V,VectorValue(0.0,0.0)) V_reg = V_φ = TestFESpace(model,reffe_scalar) U_reg = TrialFESpace(V_reg) From 480a8f4a75751ebf4694127631c26e20b6aeff24 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Tue, 12 Mar 2024 16:11:50 +1000 Subject: [PATCH 05/14] bug --- scripts/MPI/3d_inverse_homenisation_ALM.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/MPI/3d_inverse_homenisation_ALM.jl b/scripts/MPI/3d_inverse_homenisation_ALM.jl index 6b78ed2e..af24bfba 100644 --- a/scripts/MPI/3d_inverse_homenisation_ALM.jl +++ b/scripts/MPI/3d_inverse_homenisation_ALM.jl @@ -44,7 +44,7 @@ function main(mesh_partition,distribute,el_size) reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order) reffe_scalar = ReferenceFE(lagrangian,Float64,order) V = TestFESpace(model,reffe;dirichlet_tags=["origin"]) - U = TrialFESpace(_V,VectorValue(0.0,0.0,0.0)) + U = TrialFESpace(V,VectorValue(0.0,0.0,0.0)) V_reg = V_φ = TestFESpace(model,reffe_scalar) U_reg = TrialFESpace(V_reg) From 1e6bc1c37e6500cc51bbefccbf8023f7ce423422 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Tue, 12 Mar 2024 16:21:39 +1000 Subject: [PATCH 06/14] tmp adjust for testing --- scripts/MPI/3d_thermal_compliance_ALM.jl | 16 +++++++-------- .../MPI/jobs/3d_thermal_compliance_ALM.pbs | 20 +++++++++++-------- 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/scripts/MPI/3d_thermal_compliance_ALM.jl b/scripts/MPI/3d_thermal_compliance_ALM.jl index 4e84a0b3..d2e4ad12 100644 --- a/scripts/MPI/3d_thermal_compliance_ALM.jl +++ b/scripts/MPI/3d_thermal_compliance_ALM.jl @@ -3,12 +3,12 @@ using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers, using GridapSolvers: NewtonSolver -global elx = parse(Int,ARGS[1]) -global ely = parse(Int,ARGS[2]) -global elz = parse(Int,ARGS[3]) -global Px = parse(Int,ARGS[4]) -global Py = parse(Int,ARGS[5]) -global Pz = parse(Int,ARGS[6]) +# global elx = parse(Int,ARGS[1]) +# global ely = parse(Int,ARGS[2]) +# global elz = parse(Int,ARGS[3]) +# global Px = parse(Int,ARGS[4]) +# global Py = parse(Int,ARGS[5]) +# global Pz = parse(Int,ARGS[6]) """ (MPI) Minimum thermal compliance with augmented Lagrangian method in 3D. @@ -119,8 +119,8 @@ function main(mesh_partition,distribute,el_size,coef,step) end with_mpi() do distribute - mesh_partition = (Px,Py,Pz) - el_size = (elx,ely,elz) + mesh_partition = (5,5,5)#(Px,Py,Pz) + el_size = (150,150,150)#(elx,ely,elz) all_solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true -ksp_converged_reason -ksp_rtol 1.0e-12" diff --git a/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs b/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs index a1424958..ef1ff771 100644 --- a/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs +++ b/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs @@ -12,11 +12,15 @@ PROJECT_DIR=$HOME/LevelSetTopOpt/ julia --project=$PROJECT_DIR -e "using Pkg; Pkg.precompile()" -mpirun --hostfile $PBS_NODEFILE julia --project=$PROJECT_DIR --check-bounds no -O3 --compiled-modules=no \ - $PROJECT_DIR/scripts/MPI/3d_thermal_compliance_ALM.jl \ - 150 \ - 150 \ - 150 \ - 5 \ - 5 \ - 5 \ No newline at end of file +mpirun --hostfile $PBS_NODEFILE \ + julia --project=$PROJECT_DIR --check-bounds no -O3 --compiled-modules=no \ + $PROJECT_DIR/scripts/MPI/3d_thermal_compliance_ALM.jl + +# mpirun --hostfile $PBS_NODEFILE julia --project=$PROJECT_DIR --check-bounds no -O3 --compiled-modules=no \ +# $PROJECT_DIR/scripts/MPI/3d_thermal_compliance_ALM.jl \ +# 150 \ +# 150 \ +# 150 \ +# 5 \ +# 5 \ +# 5 \ No newline at end of file From 4aed4ed5f26aa73c8344cd99e80298288dfff8c3 Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Wed, 13 Mar 2024 09:59:30 +1000 Subject: [PATCH 07/14] fixes --- scripts/MPI/3d_elastic_compliance_ALM.jl | 1 - scripts/MPI/3d_inverter_ALM.jl | 5 ++--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/scripts/MPI/3d_elastic_compliance_ALM.jl b/scripts/MPI/3d_elastic_compliance_ALM.jl index 1a706773..3e5f6f24 100644 --- a/scripts/MPI/3d_elastic_compliance_ALM.jl +++ b/scripts/MPI/3d_elastic_compliance_ALM.jl @@ -64,7 +64,6 @@ function main(mesh_partition,distribute,el_size) a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N - res(u,v,φ,dΩ,dΓ_N) = a(u,v,φ,dΩ,dΓ_N) - l(v,φ,dΩ,dΓ_N) ## Optimisation functionals J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ diff --git a/scripts/MPI/3d_inverter_ALM.jl b/scripts/MPI/3d_inverter_ALM.jl index 7f5e2534..4fff7765 100644 --- a/scripts/MPI/3d_inverter_ALM.jl +++ b/scripts/MPI/3d_inverter_ALM.jl @@ -117,9 +117,8 @@ function main(mesh_partition,distribute,el_size) ) ## Optimiser - ls_enabled=false # Setting to true will use a line search instead of oscillation detection - optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,α_min=0.4,ls_enabled, - verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out]) + optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out]) for (it, uh, φh) in optimiser write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) write_history(path*"/history.txt",optimiser.history;ranks=ranks) From 7f8003a238675ba42edb085e2aa15dda2fd9fe93 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 14 Mar 2024 09:35:59 +1100 Subject: [PATCH 08/14] Bugfix: RepeatingAffineFEMap returned wrong U-assembler --- src/ChainRules.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/ChainRules.jl b/src/ChainRules.jl index 900d35db..a3ba34f6 100644 --- a/src/ChainRules.jl +++ b/src/ChainRules.jl @@ -705,10 +705,15 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap @check nblocks == length(l) spaces_0 = (U0,V0) + assem_U0 = assem_U + biform = IntegrandWithMeasure(a,dΩ) liforms = map(li -> IntegrandWithMeasure(li,dΩ),l) U, V = repeat_spaces(nblocks,U0,V0) spaces = (U,V,V_φ,U_reg) + assem_U = SparseMatrixAssembler( + get_matrix_type(assem_U0),get_vector_type(assem_V0),U,V,get_assembly_strategy(assem_U0) + ) ## Pullback cache uhd = zero(U0) @@ -802,9 +807,9 @@ function update_adjoint_caches!(φ_to_u::RepeatingAffineFEStateMap,uh,φh) return φ_to_u.adj_caches end -function adjoint_solve!(φ_to_u::RepeatingAffineFEStateMap,du::AbstractVector) +function adjoint_solve!(φ_to_u::RepeatingAffineFEStateMap,du::AbstractBlockVector) adjoint_ns, _, adjoint_x, _ = φ_to_u.adj_caches - map(blocks(adjoint_x),du) do xi, dui + map(blocks(adjoint_x),blocks(du)) do xi, dui solve!(xi,adjoint_ns,dui) end return adjoint_x From 30837bf457cd1e23d1c6079b3185c8044e12d454 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 14 Mar 2024 09:47:37 +1100 Subject: [PATCH 09/14] Minor --- scripts/Serial/inverse_homenisation_ALM.jl | 2 +- src/ChainRules.jl | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl index c480820d..f6b5ccbf 100644 --- a/scripts/Serial/inverse_homenisation_ALM.jl +++ b/scripts/Serial/inverse_homenisation_ALM.jl @@ -25,7 +25,7 @@ function main() α_coeff = 4 vf = 0.5 path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM" - mkdir(path) + !isdir(path) && mkdir(path) ## FE Setup model = CartesianDiscreteModel(dom,el_size,isperiodic=(true,true)) diff --git a/src/ChainRules.jl b/src/ChainRules.jl index a3ba34f6..7a6b0c3f 100644 --- a/src/ChainRules.jl +++ b/src/ChainRules.jl @@ -712,7 +712,7 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap U, V = repeat_spaces(nblocks,U0,V0) spaces = (U,V,V_φ,U_reg) assem_U = SparseMatrixAssembler( - get_matrix_type(assem_U0),get_vector_type(assem_V0),U,V,get_assembly_strategy(assem_U0) + get_matrix_type(assem_U0),get_vector_type(assem_U0),U,V,FESpaces.get_assembly_strategy(assem_U0) ) ## Pullback cache @@ -725,12 +725,12 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap plb_caches = (dudφ_vec,assem_deriv) ## Forward cache - K = assemble_matrix((u,v) -> biform(u,v,φh),assem_U,U0,V0) + K = assemble_matrix((u,v) -> biform(u,v,φh),assem_U0,U0,V0) b = allocate_in_range(K); fill!(b,zero(eltype(b))) b0 = allocate_in_range(K); fill!(b0,zero(eltype(b0))) x = mortar(map(i -> allocate_in_domain(K), 1:nblocks)); fill!(x,zero(eltype(x))) ns = numerical_setup(symbolic_setup(ls,K),K) - fwd_caches = (ns,K,b,x,uhd,assem_U,b0) + fwd_caches = (ns,K,b,x,uhd,assem_U0,b0,assem_U) ## Adjoint cache adjoint_K = assemble_matrix((u,v)->biform(v,u,φh),assem_adjoint,V0,U0) @@ -764,26 +764,26 @@ end get_state(m::RepeatingAffineFEStateMap) = FEFunction(get_trial_space(m),m.fwd_caches[4]) get_measure(m::RepeatingAffineFEStateMap) = m.biform.dΩ get_spaces(m::RepeatingAffineFEStateMap) = m.spaces -get_assemblers(m::RepeatingAffineFEStateMap) = (m.fwd_caches[6],m.plb_caches[2],m.adj_caches[4]) +get_assemblers(m::RepeatingAffineFEStateMap) = (m.fwd_caches[8],m.plb_caches[2],m.adj_caches[4]) 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_U, 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_U,U0,V0) + assemble_matrix!(a_fwd,K,assem_U0,U0,V0) numerical_setup!(ns,K) l0_fwd(v) = a_fwd(uhd,v) - assemble_vector!(l0_fwd,b0,assem_U,V0) + assemble_vector!(l0_fwd,b0,assem_U0,V0) rmul!(b0,-1) v = get_fe_basis(V0) map(blocks(x),liforms) do xi, li copy!(b,b0) vecdata = collect_cell_vector(V0,li(v,φh)) - assemble_vector_add!(b,assem_U,vecdata) + assemble_vector_add!(b,assem_U0,vecdata) solve!(xi,ns,b) end return x From dac40c66a5a322e2b93236b6631b8ca5f4d7d1ae Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Thu, 14 Mar 2024 13:19:59 +1000 Subject: [PATCH 10/14] clean up and fix --- scripts/MPI/3d_inverse_homenisation_ALM.jl | 18 +++++++----------- scripts/Serial/inverse_homenisation_ALM.jl | 14 ++++++-------- src/ChainRules.jl | 4 +++- 3 files changed, 16 insertions(+), 20 deletions(-) diff --git a/scripts/MPI/3d_inverse_homenisation_ALM.jl b/scripts/MPI/3d_inverse_homenisation_ALM.jl index af24bfba..b57d1785 100644 --- a/scripts/MPI/3d_inverse_homenisation_ALM.jl +++ b/scripts/MPI/3d_inverse_homenisation_ALM.jl @@ -67,16 +67,12 @@ function main(mesh_partition,distribute,el_size) l = [(v,φ,dΩ) -> ∫(-(I ∘ φ) * C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:6] ## Optimisation functionals - _C(C,ε_p,ε_q) = C ⊙ ε_p ⊙ ε_q; - - _K(C,u,εᴹ) = (_C(C,ε(u[1])+εᴹ[1],εᴹ[1]) + _C(C,ε(u[2])+εᴹ[2],εᴹ[2]) + _C(C,ε(u[3])+εᴹ[3],εᴹ[3]) + - 2(_C(C,ε(u[1])+εᴹ[1],εᴹ[2]) + _C(C,ε(u[1])+εᴹ[1],εᴹ[3]) + _C(C,ε(u[2])+εᴹ[2],εᴹ[3])))/9 - - _v_K(C,u,εᴹ) = (_C(C,ε(u[1])+εᴹ[1],ε(u[1])+εᴹ[1]) + _C(C,ε(u[2])+εᴹ[2],ε(u[2])+εᴹ[2]) + _C(C,ε(u[3])+εᴹ[3],ε(u[3])+εᴹ[3]) + - 2(_C(C,ε(u[1])+εᴹ[1],ε(u[2])+εᴹ[2]) + _C(C,ε(u[1])+εᴹ[1],ε(u[3])+εᴹ[3]) + _C(C,ε(u[2])+εᴹ[2],ε(u[3])+εᴹ[3])))/9 - - J(u,φ,dΩ) = ∫(-(I ∘ φ)*_K(C,u,εᴹ))dΩ - dJ(q,u,φ,dΩ) = ∫(_v_K(C,u,εᴹ)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ,dΩ) = ∫(q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ,dΩ) = -1/9*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+Cᴴ(3,3,u,φ,dΩ)+ + 2*(Cᴴ(1,2,u,φ,dΩ)+Cᴴ(1,3,u,φ,dΩ)+Cᴴ(2,3,u,φ,dΩ))) + dκ(q,u,φ,dΩ) = 1/9*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+dCᴴ(3,3,q,u,φ,dΩ)+ + 2*(dCᴴ(1,2,q,u,φ,dΩ)+dCᴴ(1,3,q,u,φ,dΩ)+dCᴴ(2,3,q,u,φ,dΩ))) Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ @@ -95,7 +91,7 @@ function main(mesh_partition,distribute,el_size) assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), ls = solver, adjoint_ls = solver ) - pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dJ=dJ,analytic_dC=[dVol]) + pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems α = α_coeff*maximum(el_Δ) diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl index c480820d..4dcf8b1f 100644 --- a/scripts/Serial/inverse_homenisation_ALM.jl +++ b/scripts/Serial/inverse_homenisation_ALM.jl @@ -24,7 +24,7 @@ function main() η_coeff = 2 α_coeff = 4 vf = 0.5 - path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM" + path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM_new" mkdir(path) ## FE Setup @@ -62,12 +62,10 @@ function main() l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3] ## Optimisation functionals - _C(C,ε_p,ε_q) = C ⊙ ε_p ⊙ ε_q - _K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],εᴹ[2]))/4 - _v_K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],ε(u1)+εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],ε(u2)+εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],ε(u2)+εᴹ[2]))/4 - - J(u,φ,dΩ) = ∫(-(I ∘ φ)*_K(C,u,εᴹ))dΩ - dJ(q,u,φ,dΩ) = ∫(_v_K(C,u,εᴹ)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ,dΩ) = ∫(q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ)) + dκ(q,u,φ,dΩ) = 1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ)) Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ @@ -76,7 +74,7 @@ function main() ## Setup solver and FE operators state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ) - pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dJ=dJ,analytic_dC=[dVol]) + pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems α = α_coeff*maximum(el_Δ) diff --git a/src/ChainRules.jl b/src/ChainRules.jl index 900d35db..d9387dbe 100644 --- a/src/ChainRules.jl +++ b/src/ChainRules.jl @@ -49,7 +49,9 @@ function Gridap.gradient(F::IntegrandWithMeasure,uh::Vector,K::Int) local_fields = map(local_views,uh) |> to_parray_of_arrays local_measures = map(local_views,F.dΩ) |> to_parray_of_arrays contribs = map(local_measures,local_fields) do dΩ,lf - _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...) + # TODO: Remove second term below, this is a fix for the problem discussed in + # https://github.com/zjwegert/LSTO_Distributed/issues/46 + _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...) #+ ∑(∫(0)dΩ[i] for i = 1:length(dΩ)) return Gridap.Fields.gradient(_f,lf[K]) end return DistributedDomainContribution(contribs) From 274a03f814381d10c90d5755aa3b6ede091e6779 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Thu, 14 Mar 2024 13:28:25 +1000 Subject: [PATCH 11/14] fix name --- scripts/Serial/inverse_homenisation_ALM.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl index 4dcf8b1f..5090f339 100644 --- a/scripts/Serial/inverse_homenisation_ALM.jl +++ b/scripts/Serial/inverse_homenisation_ALM.jl @@ -24,7 +24,7 @@ function main() η_coeff = 2 α_coeff = 4 vf = 0.5 - path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM_new" + path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM" mkdir(path) ## FE Setup From f548e1e5d939ae7c4f7857d27689a769fbc9564f Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Thu, 14 Mar 2024 14:15:02 +1000 Subject: [PATCH 12/14] fix and add scripts for new issue --- scripts/MPI/3d_inverse_homenisation_ALM.jl | 4 +- scripts/Serial/inverse_homenisation_ALM.jl | 4 +- scripts/_dev/bug_issue_46/thermal_MPI.jl | 91 +++++++++++++++++++++ scripts/_dev/bug_issue_46/thermal_serial.jl | 87 ++++++++++++++++++++ 4 files changed, 182 insertions(+), 4 deletions(-) create mode 100644 scripts/_dev/bug_issue_46/thermal_MPI.jl create mode 100644 scripts/_dev/bug_issue_46/thermal_serial.jl diff --git a/scripts/MPI/3d_inverse_homenisation_ALM.jl b/scripts/MPI/3d_inverse_homenisation_ALM.jl index b57d1785..488ca57e 100644 --- a/scripts/MPI/3d_inverse_homenisation_ALM.jl +++ b/scripts/MPI/3d_inverse_homenisation_ALM.jl @@ -68,10 +68,10 @@ function main(mesh_partition,distribute,el_size) ## Optimisation functionals Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ - dCᴴ(r,s,q,u,φ,dΩ) = ∫(q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ κ(u,φ,dΩ) = -1/9*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+Cᴴ(3,3,u,φ,dΩ)+ 2*(Cᴴ(1,2,u,φ,dΩ)+Cᴴ(1,3,u,φ,dΩ)+Cᴴ(2,3,u,φ,dΩ))) - dκ(q,u,φ,dΩ) = 1/9*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+dCᴴ(3,3,q,u,φ,dΩ)+ + dκ(q,u,φ,dΩ) = -1/9*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+dCᴴ(3,3,q,u,φ,dΩ)+ 2*(dCᴴ(1,2,q,u,φ,dΩ)+dCᴴ(1,3,q,u,φ,dΩ)+dCᴴ(2,3,q,u,φ,dΩ))) Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl index 5090f339..d4d60aa8 100644 --- a/scripts/Serial/inverse_homenisation_ALM.jl +++ b/scripts/Serial/inverse_homenisation_ALM.jl @@ -63,9 +63,9 @@ function main() ## Optimisation functionals Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ - dCᴴ(r,s,q,u,φ,dΩ) = ∫(q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ)) - dκ(q,u,φ,dΩ) = 1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ)) + dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ)) Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ diff --git a/scripts/_dev/bug_issue_46/thermal_MPI.jl b/scripts/_dev/bug_issue_46/thermal_MPI.jl new file mode 100644 index 00000000..2e4596ee --- /dev/null +++ b/scripts/_dev/bug_issue_46/thermal_MPI.jl @@ -0,0 +1,91 @@ +using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays, + SparseMatricesCSR + +function main(mesh_partition,distribute,use_l::Bool) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + # FE parameters + order = 1 # Finite element order + xmax = ymax = 1.0 # Domain size + dom = (0,xmax,0,ymax) # Bounding domain + el_size = (200,200) # Mesh partition size + prop_Γ_N = 0.2 # Γ_N size parameter + prop_Γ_D = 0.2 # Γ_D size parameter + f_Γ_N(x) = (x[1] ≈ xmax && # Γ_N indicator function + ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + eps()) + f_Γ_D(x) = (x[1] ≈ 0.0 && # Γ_D indicator function + (x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps())) # @$\lvert\lvert$@ + # FD parameters + γ = 0.1 # HJ eqn time step coeff + γ_reinit = 0.5 # Reinit. eqn time step coeff + max_steps = floor(Int,minimum(el_size)/10) # Max steps for advection + tol = 1/(5order^2)/minimum(el_size) # Advection tolerance + # Problem parameters + κ = 1 # Diffusivity + g = 1 # Heat flow in + vf = 0.4 # Volume fraction constraint + lsf_func = initial_lsf(4,0.2) # Initial level set function + iter_mod = 10 # VTK Output modulo + path = "./results/AD_tut1_MPI_lobj_$use_l/" # Output path + i_am_main(ranks) && mkpath(path) # Create path + # Model + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); + update_labels!(1,model,f_Γ_D,"Gamma_D") + update_labels!(2,model,f_Γ_N,"Gamma_N") + # Triangulation and measures + Ω = Triangulation(model) + Γ_N = BoundaryTriangulation(model,tags="Gamma_N") + dΩ = Measure(Ω,2*order) + dΓ_N = Measure(Γ_N,2*order) + # Spaces + reffe = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,0.0) + V_φ = TestFESpace(model,reffe) + V_reg = TestFESpace(model,reffe;dirichlet_tags=["Gamma_N"]) + U_reg = TrialFESpace(V_reg,0) + # Level set and interpolator + φh = interpolate(lsf_func,V_φ) + interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model))) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + # Weak formulation + a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(g*v)dΓ_N + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + # Objective and constraints + J(u,φ,dΩ,dΓ_N) = if use_l + ∫(g*u)dΓ_N #+ ∫(0)dΩ + else + ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + end + dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + vol_D = sum(∫(1)dΩ) + C(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + pcfs = PDEConstrainedFunctionals(J,[C],state_map) + # Velocity extension + α = 4*maximum(get_el_Δ(model)) + a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ + vel_ext = VelocityExtension(a_hilb, U_reg, V_reg) + # Finite difference scheme + scheme = FirstOrderStencil(length(el_size),Float64) + stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + # Optimiser + optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) + # Solve + for (it,uh,φh) in optimiser + data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] + iszero(it % iter_mod) && (writevtk(Ω,path*"out$it",cellfields=data);GC.gc()) + write_history(path*"/history.txt",get_history(optimiser);ranks) + end + # Final structure + it = get_history(optimiser).niter; uh = get_state(pcfs) + writevtk(Ω,path*"out$it",cellfields=["φ"=>φh, + "H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +end + +with_mpi() do distribute + mesh_partition = (2,2) + main(mesh_partition,distribute,false) + main(mesh_partition,distribute,true) # <- this should crash +end \ No newline at end of file diff --git a/scripts/_dev/bug_issue_46/thermal_serial.jl b/scripts/_dev/bug_issue_46/thermal_serial.jl new file mode 100644 index 00000000..e9c29a6c --- /dev/null +++ b/scripts/_dev/bug_issue_46/thermal_serial.jl @@ -0,0 +1,87 @@ +using LevelSetTopOpt, Gridap + +function main(use_l::Bool) + # FE parameters + order = 1 # Finite element order + xmax = ymax = 1.0 # Domain size + dom = (0,xmax,0,ymax) # Bounding domain + el_size = (200,200) # Mesh partition size + prop_Γ_N = 0.2 # Γ_N size parameter + prop_Γ_D = 0.2 # Γ_D size parameter + f_Γ_N(x) = (x[1] ≈ xmax && # Γ_N indicator function + ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + eps()) + f_Γ_D(x) = (x[1] ≈ 0.0 && # Γ_D indicator function + (x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps())) + # FD parameters + γ = 0.1 # HJ eqn time step coeff + γ_reinit = 0.5 # Reinit. eqn time step coeff + max_steps = floor(Int,minimum(el_size)/10) # Max steps for advection + tol = 1/(5order^2)/minimum(el_size) # Reinitialisation tolerance + # Problem parameters + κ = 1 # Diffusivity + g = 1 # Heat flow in + vf = 0.4 # Volume fraction constraint + lsf_func = initial_lsf(4,0.2) # Initial level set function + iter_mod = 10 # VTK Output modulo + path = "./results/tut1_lobj_$use_l/" # Output path + mkpath(path) # Create path + # Model + model = CartesianDiscreteModel(dom,el_size); + update_labels!(1,model,f_Γ_D,"Gamma_D") + update_labels!(2,model,f_Γ_N,"Gamma_N") + # Triangulation and measures + Ω = Triangulation(model) + Γ_N = BoundaryTriangulation(model,tags="Gamma_N") + dΩ = Measure(Ω,2*order) + dΓ_N = Measure(Γ_N,2*order) + # Spaces + reffe = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,0.0) + V_φ = TestFESpace(model,reffe) + V_reg = TestFESpace(model,reffe;dirichlet_tags=["Gamma_N"]) + U_reg = TrialFESpace(V_reg,0) + # Level set and interpolator + φh = interpolate(lsf_func,V_φ) + interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model))) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + # Weak formulation + a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(g*v)dΓ_N + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + # Objective and constraints + J(u,φ,dΩ,dΓ_N) = if use_l + ∫(g*u)dΓ_N + else + ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + end + dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + vol_D = sum(∫(1)dΩ) + C(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + pcfs = PDEConstrainedFunctionals(J,[C],state_map, + analytic_dJ=dJ,analytic_dC=[dC]) + # Velocity extension + α = 4*maximum(get_el_Δ(model)) + a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ + vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + # Finite difference scheme + scheme = FirstOrderStencil(length(el_size),Float64) + stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps) + # Optimiser + optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; + γ,γ_reinit,verbose=true,constraint_names=[:Vol]) + # Solve + for (it,uh,φh) in optimiser + data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh] + iszero(it % iter_mod) && (writevtk(Ω,path*"out$it",cellfields=data);GC.gc()) + write_history(path*"/history.txt",get_history(optimiser)) + end + # Final structure + it = get_history(optimiser).niter; uh = get_state(pcfs) + writevtk(Ω,path*"out$it",cellfields=["φ"=>φh, + "H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +end + +main(false) +main(true) \ No newline at end of file From 67922d71f0319e5c673fee425357003f3d4ab88e Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Thu, 14 Mar 2024 14:33:31 +1000 Subject: [PATCH 13/14] Added new benchmark --- scripts/Benchmarks/benchmark.jl | 2 +- src/Benchmarks.jl | 19 +++++++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/scripts/Benchmarks/benchmark.jl b/scripts/Benchmarks/benchmark.jl index d91cbe6f..9abd59ee 100644 --- a/scripts/Benchmarks/benchmark.jl +++ b/scripts/Benchmarks/benchmark.jl @@ -439,7 +439,7 @@ with_mpi() do distribute bopt0 = i_am_main(ranks) && zeros(NREPS) end if occursin("bopt1",BMARK_TYPE) - bopt1 = benchmark_optimizer(optim, 1, ranks; nreps=NREPS) + bopt1 = benchmark_single_iteration(optim, ranks; nreps = NREPS) else bopt1 = i_am_main(ranks) && zeros(NREPS) end diff --git a/src/Benchmarks.jl b/src/Benchmarks.jl index d59ff5bd..87cc4a6e 100644 --- a/src/Benchmarks.jl +++ b/src/Benchmarks.jl @@ -79,6 +79,25 @@ function benchmark_optimizer(m::Optimiser, niter, ranks; nreps = 10) return benchmark(f, (m,), ranks; nreps, reset! = opt_reset!) end +""" + benchmark_single_iteration(m::Optimiser, ranks; nreps) + +Given an optimiser `m`, benchmark a single iteration after 0th iteration. +""" +function benchmark_single_iteration(m::Optimiser, ranks; nreps = 10) + _, state = iterate(m) + function f(m) + iterate(m, state) + end + + φ0 = copy(get_free_dof_values(m.φ0)) + function opt_reset!(m::Optimiser) + copy!(get_free_dof_values(m.φ0), φ0) + reset!(get_history(m)) + end + return benchmark(f, (m,), ranks; nreps, reset! = opt_reset!) +end + """ benchmark_forward_problem(m::AbstractFEStateMap, φh, ranks; nreps) From 9508a1c6da1d13bae1bc22a73590ee41c8f7bf1c Mon Sep 17 00:00:00 2001 From: "Z J Wegert (Workstation)" Date: Mon, 18 Mar 2024 14:36:13 +1000 Subject: [PATCH 14/14] Testing --- scripts/MPI/3d_thermal_compliance_ALM.jl | 12 ++++++------ scripts/MPI/inverse_homenisation_ALM.jl | 12 +++++------- scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs | 2 +- 3 files changed, 12 insertions(+), 14 deletions(-) diff --git a/scripts/MPI/3d_thermal_compliance_ALM.jl b/scripts/MPI/3d_thermal_compliance_ALM.jl index d2e4ad12..cae68f2d 100644 --- a/scripts/MPI/3d_thermal_compliance_ALM.jl +++ b/scripts/MPI/3d_thermal_compliance_ALM.jl @@ -73,7 +73,7 @@ function main(mesh_partition,distribute,el_size,coef,step) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N + l(v,φ,dΩ,dΓ_N) = ∫(10v)dΓ_N ## Optimisation functionals J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ @@ -126,12 +126,12 @@ with_mpi() do distribute GridapPETSc.with(args=split(all_solver_options)) do main(mesh_partition,distribute,el_size,5,10) - main(mesh_partition,distribute,el_size,2,10) + # main(mesh_partition,distribute,el_size,2,10) - main(mesh_partition,distribute,el_size,5,5) - main(mesh_partition,distribute,el_size,2,5) + # main(mesh_partition,distribute,el_size,5,5) + # main(mesh_partition,distribute,el_size,2,5) - main(mesh_partition,distribute,el_size,5,3) - main(mesh_partition,distribute,el_size,2,3) + # main(mesh_partition,distribute,el_size,5,3) + # main(mesh_partition,distribute,el_size,2,3) end end; \ No newline at end of file diff --git a/scripts/MPI/inverse_homenisation_ALM.jl b/scripts/MPI/inverse_homenisation_ALM.jl index ed9c279f..e5399ae1 100644 --- a/scripts/MPI/inverse_homenisation_ALM.jl +++ b/scripts/MPI/inverse_homenisation_ALM.jl @@ -64,12 +64,10 @@ function main(mesh_partition,distribute,el_size) l = [(v,φ,dΩ) -> ∫(-(I ∘ φ) * C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3] ## Optimisation functionals - _C(C,ε_p,ε_q) = C ⊙ ε_p ⊙ ε_q; - _K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],εᴹ[2]))/4 - _v_K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],ε(u1)+εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],ε(u2)+εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],ε(u2)+εᴹ[2]))/4 - - J(u,φ,dΩ) = ∫(-(I ∘ φ)*_K(C,u,εᴹ))dΩ - dJ(q,u,φ,dΩ) = ∫(_v_K(C,u,εᴹ)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ)) + dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ)) Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ @@ -88,7 +86,7 @@ function main(mesh_partition,distribute,el_size) assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), ls = solver,adjoint_ls = solver ) - pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dJ=dJ,analytic_dC=[dVol]) + pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems α = α_coeff*maximum(el_Δ) diff --git a/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs b/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs index ef1ff771..d71b2b8f 100644 --- a/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs +++ b/scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs @@ -4,7 +4,7 @@ #PBS -N 3d_thermal_compliance_ALM #PBS -l select=125:ncpus=1:mpiprocs=1:ompthreads=1:mem=8GB:cputype=6140 -#PBS -l walltime=200:00:00 +#PBS -l walltime=23:00:00 #PBS -j oe source $HOME/hpc-environments-main/lyra/load-ompi.sh