diff --git a/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl b/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl index 0ddd1da4..cc428534 100644 --- a/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl +++ b/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl @@ -61,9 +61,9 @@ function main(mesh_partition,distribute,el_size) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ ## Material properties - _E = 1000 + E = 1000 ν = 0.3 - μ, λ = _E/(2*(1 + ν)), _E*ν/((1 + ν)*(1 - 2*ν)) + μ, λ = E/(2*(1 + ν)), E*ν/((1 + ν)*(1 - 2*ν)) g = VectorValue(0,0,-100) ## Neohookean hyperelastic material diff --git a/scripts/_dev/Added_for_paper/hyperelastic_compliance_neohook_ALM_COMPARE.jl b/scripts/_dev/Added_for_paper/hyperelastic_compliance_neohook_ALM_COMPARE.jl deleted file mode 100644 index 673ae38c..00000000 --- a/scripts/_dev/Added_for_paper/hyperelastic_compliance_neohook_ALM_COMPARE.jl +++ /dev/null @@ -1,193 +0,0 @@ -using Gridap, LevelSetTopOpt - -""" - (Serial) Minimum hyperelastic (neohookean) compliance with augmented Lagrangian method in 2D. -""" -function main() - ## Parameters - order = 1 - xmax,ymax=2.0,1.0 - prop_Γ_N = 0.2 - dom = (0,xmax,0,ymax) - el_size = (400,200) - γ = 0.1 - γ_reinit = 0.5 - max_steps = floor(Int,minimum(el_size)/10) - tol = 1/(5order^2)/minimum(el_size) - η_coeff = 2 - α_coeff = 4 - vf = 0.4 - path = dirname(dirname(@__DIR__))*"/results/highres_hyperelastic_compliance_neohook_ALM" - mkdir(path) - - ## FE Setup - model = CartesianDiscreteModel(dom,el_size); - el_Δ = get_el_Δ(model) - f_Γ_D(x) = (x[1] ≈ 0.0) - f_Γ_N(x) = (x[1] ≈ xmax && ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + eps()) - update_labels!(1,model,f_Γ_D,"Gamma_D") - update_labels!(2,model,f_Γ_N,"Gamma_N") - - ## Triangulations and measures - Ω = Triangulation(model) - Γ_N = BoundaryTriangulation(model,tags="Gamma_N") - dΩ = Measure(Ω,2*order) - dΓ_N = Measure(Γ_N,2*order) - vol_D = sum(∫(1)dΩ) - - ## Spaces - reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) - reffe_scalar = ReferenceFE(lagrangian,Float64,order) - V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) - U = TrialFESpace(V,VectorValue(0.0,0.0)) - V_φ = TestFESpace(model,reffe_scalar) - V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) - U_reg = TrialFESpace(V_reg,0) - - ## Create FE functions - φh = interpolate(initial_lsf(4,0.2),V_φ) - - ## Interpolation and weak form - interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) - I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ - - ## Material properties - _E = 1000 - ν = 0.3 - μ, λ = _E/(2*(1 + ν)), _E*ν/((1 + ν)*(1 - 2*ν)) - g = VectorValue(0,-10) - - ## Neohookean hyperelastic material - # Deformation gradient - F(∇u) = one(∇u) + ∇u' - J(F) = sqrt(det(C(F))) - - # Derivative of green Strain - dE(∇du,∇u) = 0.5*( ∇du⋅F(∇u) + (∇du⋅F(∇u))' ) - - # Right Caughy-green deformation tensor - C(F) = (F')⋅F - - # Constitutive law (Neo hookean) - function S(∇u) - Cinv = inv(C(F(∇u))) - μ*(one(∇u)-Cinv) + λ*log(J(F(∇u)))*Cinv - end - - # Cauchy stress tensor and residual - σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))' - res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N - - ## Optimisation functionals - Obj(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ - Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; - dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ - - ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) - - ## Setup solver and FE operators - state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) - pcfs = PDEConstrainedFunctionals(Obj,[Vol],state_map,analytic_dC=[dVol]) - - ## Hilbertian extension-regularisation problems - α = α_coeff*maximum(el_Δ) - a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ - vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) - - ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; - γ,γ_reinit,verbose=true,constraint_names=[:Vol]) - 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) - end - it = get_history(optimiser).niter; uh = get_state(pcfs) - write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1) -end - -main() - -function main_LE_comparison() - ## Parameters - order = 1 - xmax,ymax=(2.0,1.0) - prop_Γ_N = 0.2 - dom = (0,xmax,0,ymax) - el_size = (400,200) - γ = 0.1 - γ_reinit = 0.5 - max_steps = floor(Int,minimum(el_size)/10) - tol = 1/(5order^2)/minimum(el_size) - C = isotropic_elast_tensor(2,1000,0.3) - η_coeff = 2 - α_coeff = 4 - vf = 0.4 - g = VectorValue(0,-10) - path = dirname(dirname(@__DIR__))*"/results/highres_LE_comparison" - mkdir(path) - - ## FE Setup - model = CartesianDiscreteModel(dom,el_size) - el_Δ = get_el_Δ(model) - f_Γ_D(x) = (x[1] ≈ 0.0) - f_Γ_N(x) = (x[1] ≈ xmax && ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + eps()) - update_labels!(1,model,f_Γ_D,"Gamma_D") - update_labels!(2,model,f_Γ_N,"Gamma_N") - - ## Triangulations and measures - Ω = Triangulation(model) - Γ_N = BoundaryTriangulation(model,tags="Gamma_N") - dΩ = Measure(Ω,2*order) - dΓ_N = Measure(Γ_N,2*order) - vol_D = sum(∫(1)dΩ) - - ## Spaces - reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) - reffe_scalar = ReferenceFE(lagrangian,Float64,order) - V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) - U = TrialFESpace(V,VectorValue(0.0,0.0)) - V_φ = TestFESpace(model,reffe_scalar) - V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) - U_reg = TrialFESpace(V_reg,0) - - ## Create FE functions - φh = interpolate(initial_lsf(4,0.2),V_φ) - - ## Interpolation and weak form - interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) - I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ - - a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N - - ## Optimisation functionals - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ - dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; - Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; - dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ - - ## Finite difference solver and level set function - stencil = AdvectionStencil(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) - - ## Setup solver and FE operators - state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) - pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) - - ## Hilbertian extension-regularisation problems - α = α_coeff*maximum(el_Δ) - a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; - vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) - - ## Optimiser - optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh; - γ,γ_reinit,verbose=true,constraint_names=[:Vol]) - 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) - end - it = get_history(optimiser).niter; uh = get_state(pcfs) - write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1) -end - -main_LE_comparison() \ No newline at end of file diff --git a/scripts/_dev/Added_for_paper/thermal_MPI.jl b/scripts/_dev/Added_for_paper/thermal_MPI.jl deleted file mode 100644 index 446c0a54..00000000 --- a/scripts/_dev/Added_for_paper/thermal_MPI.jl +++ /dev/null @@ -1,88 +0,0 @@ -using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays, - SparseMatricesCSR - -function main(mesh_partition,distribute) - 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/tut1_MPI_FIXEDSign/" # 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 - solver = PETScLinearSolver() - 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=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) -end \ No newline at end of file diff --git a/scripts/_dev/Added_for_paper/thermal_PETSc.jl b/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl similarity index 50% rename from scripts/_dev/Added_for_paper/thermal_PETSc.jl rename to scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl index 9b1383cf..b3bd339d 100644 --- a/scripts/_dev/Added_for_paper/thermal_PETSc.jl +++ b/scripts/_dev/Paper_scripts/elast_comp_MPI_3D.jl @@ -1,99 +1,104 @@ -using LevelSetTopOpt, Gridap, GridapPETSc, SparseMatricesCSR - -function main() - # 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_PETSc/" # 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 - Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} - Tv = Vector{PetscScalar} - solver = PETScLinearSolver() - state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; - assem_U = SparseMatrixAssembler(Tm,Tv,U,V), - assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), - assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), - ls = solver,adjoint_ls = solver - ) - # 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; - assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg), - ls = solver - ) - # 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 - -solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true - -ksp_converged_reason -ksp_rtol 1.0e-12" -GridapPETSc.with(args=split(solver_options)) do - main() +using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays, + SparseMatricesCSR + +function main(mesh_partition,distribute) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + # FE parameters + order = 1 # Finite element order + xmax,ymax,zmax = (2.0,1.0,1.0) # Domain size + dom = (0,xmax,0,ymax,0,zmax) # Bounding domain + el_size = (160,80,80) # Mesh partition size + prop_Γ_N = 0.2 # Γ_N 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())&& + (zmax/2-zmax*prop_Γ_N/2 - eps() <= x[3] <= zmax/2+zmax*prop_Γ_N/2 + eps()) + f_Γ_D(x) = (x[1] ≈ 0.0) # Γ_D indicator function + # 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 + C = isotropic_elast_tensor(3,1.,0.3) # Stiffness tensor + g = VectorValue(0,0,-1) # Applied load on Γ_N + 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/elast_comp_MPI_3D/" # 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,VectorValue{3,Float64},order) + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,VectorValue(0.0,0.0,0.0)) + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar;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 ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N + Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} + Tv = Vector{PetscScalar} + solver = ElasticitySolver(V) + state_map = AffineFEStateMap( + a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + assem_U = SparseMatrixAssembler(Tm,Tv,U,V), + assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), + assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), + ls = solver,adjoint_ls = solver + ) + # Objective and constraints + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ + dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + vol_D = sum(∫(1)dΩ) + C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + pcfs = PDEConstrainedFunctionals(J,[C1],state_map, + analytic_dJ=dJ,analytic_dC=[dC1]) + # 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; + assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg), + ls = PETScLinearSolver() + ) + # 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)) + # 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) + 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),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +end + +with_mpi() do distribute + mesh_partition = (5,5,5) + solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true + -ksp_converged_reason -ksp_rtol 1.0e-12" + GridapPETSc.with(args=split(solver_options)) do + main(mesh_partition,distribute) + end end \ No newline at end of file diff --git a/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl b/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl new file mode 100644 index 00000000..919e046e --- /dev/null +++ b/scripts/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl @@ -0,0 +1,118 @@ +using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays, + SparseMatricesCSR + +function main(mesh_partition,distribute) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + # FE parameters + order = 1 # Finite element order + xmax,ymax,zmax = (2.0,1.0,1.0) # Domain size + dom = (0,xmax,0,ymax,0,zmax) # Bounding domain + el_size = (160,80,80) # Mesh partition size + prop_Γ_N = 0.2 # Γ_N 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())&& + (zmax/2-zmax*prop_Γ_N/2 - eps() <= x[3] <= zmax/2+zmax*prop_Γ_N/2 + eps()) + f_Γ_D(x) = (x[1] ≈ 0.0) # Γ_D indicator function + # 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 + E = 1000 # Young's modulus + ν = 0.3 # Poisson's ratio + μ, λ = E/(2*(1 + ν)), E*ν/((1 + ν)*(1 - 2*ν)) # Lame constants + g = VectorValue(0,0,-1) # Applied load on Γ_N + 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/hyperelast_comp_MPI_3D/" # 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,VectorValue{3,Float64},order) + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,VectorValue(0.0,0.0,0.0)) + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar;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 + ## Piola-Kirchhoff tensor + function S(∇u) + Cinv = inv(C(F(∇u))) + μ*(one(∇u)-Cinv) + λ*log(J(F(∇u)))*Cinv + end + ## Derivative of Green strain + dE(∇du,∇u) = 0.5*( ∇du⋅F(∇u) + (∇du⋅F(∇u))' ) + ## Right Cauchy-Green deformation tensor + C(F) = (F')⋅F + ## Deformation gradient tensor + F(∇u) = one(∇u) + ∇u' + ## Volume change + J(F) = sqrt(det(C(F))) + ## Residual + res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE ∘ (∇(v),∇(u))) ⊙ (S ∘ ∇(u))) )*dΩ - ∫(g⋅v)dΓ_N + Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} + Tv = Vector{PetscScalar} + lin_solver = ElasticitySolver(V) + nl_solver = NewtonSolver(lin_solver;maxiter=50,rtol=10^-8,verbose=i_am_main(ranks)) + state_map = NonlinearFEStateMap( + res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + assem_U = SparseMatrixAssembler(Tm,Tv,U,V), + assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), + assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), + nls = nl_solver, adjoint_ls = lin_solver + ) + # Objective and constraints + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE ∘ (∇(u),∇(u))) ⊙ (S ∘ ∇(u))))dΩ + vol_D = sum(∫(1)dΩ) + C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + pcfs = PDEConstrainedFunctionals(J,[C1],state_map,analytic_dC=[dC1]) + # 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; + assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg), + ls = PETScLinearSolver() + ) + # 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)) + # 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) + 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),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +end + +with_mpi() do distribute + mesh_partition = (5,5,5) + solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true + -ksp_converged_reason -ksp_rtol 1.0e-12" + GridapPETSc.with(args=split(solver_options)) do + main(mesh_partition,distribute) + end +end \ No newline at end of file diff --git a/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl b/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl new file mode 100644 index 00000000..eab453d0 --- /dev/null +++ b/scripts/_dev/Paper_scripts/inverse_hom_MPI_3D.jl @@ -0,0 +1,108 @@ +using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays, + SparseMatricesCSR + +function main(mesh_partition,distribute) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + # FE parameters + order = 1 # Finite element order + xmax,ymax,zmax = (1.0,1.0,1.0) # Domain size + dom = (0,xmax,0,ymax,0,zmax) # Bounding domain + el_size = (100,100,100) # Mesh partition size + f_Γ_D(x) = iszero(x) # Γ_D indicator function + # 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 + C = isotropic_elast_tensor(3,1.,0.3) # Stiffness tensor + g = VectorValue(0,0,-1) # Applied load on Γ_N + vf = 0.4 # Volume fraction constraint + lsf_func(x) = cos(2π*x[1]) + cos(2π*x[2]) + # Initial level set function + cos(2π*x[3]) + iter_mod = 10 # VTK Output modulo + path = "./results/inverse_hom_MPI_3D/" # Output path + i_am_main(ranks) && mkpath(path) # Create path + # Model + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size,isperiodic=(true,true,true)); + update_labels!(1,model,f_Γ_D,"origin") + # Triangulation and measures + Ω = Triangulation(model) + dΩ = Measure(Ω,2*order) + # Spaces + 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)) + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar) + U_reg = TrialFESpace(V_reg) + # 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 + εᴹ = (TensorValue(1.,0.,0.,0.,0.,0.,0.,0.,0.), + TensorValue(0.,0.,0.,0.,1.,0.,0.,0.,0.), + TensorValue(0.,0.,0.,0.,0.,0.,0.,0.,1.), + TensorValue(0.,0.,0.,0.,0.,1/2,0.,1/2,0.), + TensorValue(0.,0.,1/2,0.,0.,0.,1/2,0.,0.), + TensorValue(0.,1/2,0.,1/2,0.,0.,0.,0.,0.)) + a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v))dΩ + l = [(v,φ,dΩ) -> ∫(-(I ∘ φ) * C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:6] + Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} + Tv = Vector{PetscScalar} + solver = ElasticitySolver(V) + state_map = RepeatingAffineFEStateMap( + 6,a,l,U,V,V_φ,U_reg,φh,dΩ; + assem_U = SparseMatrixAssembler(Tm,Tv,U,V), + assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), + assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), + ls = solver, adjoint_ls = solver + ) + # Objective and constraints + 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_D = sum(∫(1)dΩ) + C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + pcfs = PDEConstrainedFunctionals(κ,[C1],state_map, + analytic_dJ=dκ,analytic_dC=[dC1]) + # 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; + assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg), + ls = PETScLinearSolver() + ) + # 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)) + # Solve + for (it,uh,φh) in optimiser + data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))] + iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data) + 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),"|∇(φ)|"=>(norm ∘ ∇(φh))]) +end + +with_mpi() do distribute + mesh_partition = (5,5,5) + solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true + -ksp_converged_reason -ksp_rtol 1.0e-12" + GridapPETSc.with(args=split(solver_options)) do + main(mesh_partition,distribute) + end +end \ No newline at end of file diff --git a/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl b/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl new file mode 100644 index 00000000..0dbac675 --- /dev/null +++ b/scripts/_dev/Paper_scripts/inverter_MPI_3D.jl @@ -0,0 +1,114 @@ +using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays, + SparseMatricesCSR + +function main(mesh_partition,distribute) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + # FE parameters + order = 1 # Finite element order + xmax,ymax,zmax = (1.0,1.0,1.0) # Domain size + dom = (0,xmax,0,ymax,0,zmax) # Bounding domain + el_size = (100,100,100) # Mesh partition size + f_Γ_in(x) = (x[1] ≈ 0.0) && # Γ_in indicator function + (0.4 - eps() <= x[2] <= 0.6 + eps()) && (0.4 - eps() <= x[3] <= 0.6 + eps()) + f_Γ_out(x) = (x[1] ≈ 1.0) && # Γ_out indicator function + (0.4 - eps() <= x[2] <= 0.6 + eps()) && (0.4 - eps() <= x[3] <= 0.6 + eps()) + f_Γ_D(x) = (x[1] ≈ 0.0) && # Γ_D indicator function + (x[2] <= 0.1 || x[2] >= 0.9) && (x[3] <= 0.1 || x[3] >= 0.9) + # 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 + C = isotropic_elast_tensor(3,1.,0.3) # Stiffness tensor + g = VectorValue(0,0,-1) # Applied load on Γ_N + vf = 0.4 # Volume fraction constraint + sphere(x,(xc,yc,zc)) = -sqrt((x[1]-xc)^2+(x[2]-yc)^2+(x[3]-zc)^2) + 0.2 + lsf_func(x) = max(initial_lsf(4,0.2)(x), # Initial level set function + sphere(x,(1,0,0)),sphere(x,(1,0,1)),sphere(x,(1,1,0)),sphere(x,(1,1,1))) + lsf_func = initial_lsf(4,0.2) # Initial level set function + iter_mod = 10 # VTK Output modulo + path = "./results/inverter_MPI_3D/" # Output path + i_am_main(ranks) && mkpath(path) # Create path + # Model + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); + update_labels!(1,model,f_Γ_in,"Gamma_in") + update_labels!(2,model,f_Γ_out,"Gamma_out") + update_labels!(4,model,f_Γ_D,"Gamma_D") + # Triangulation and measures + Ω = Triangulation(model) + Γ_in = BoundaryTriangulation(model,tags="Gamma_in") + Γ_out = BoundaryTriangulation(model,tags="Gamma_out") + dΩ = Measure(Ω,2*order) + dΓ_in = Measure(Γ_in,2*order) + dΓ_out = Measure(Γ_out,2*order) + # Spaces + reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order) + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,VectorValue(0.0,0.0,0.0)) + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_in","Gamma_out"]) + U_reg = TrialFESpace(V_reg,0) + # 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Γ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in + Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} + Tv = Vector{PetscScalar} + solver = ElasticitySolver(V) + state_map = AffineFEStateMap( + a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + assem_U = SparseMatrixAssembler(Tm,Tv,U,V), + assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), + assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), + ls = solver,adjoint_ls = solver + ) + # Objective and constraints + e₁ = VectorValue(1,0,0) + vol_Γ_in = sum(∫(1)dΓ_in) + vol_Γ_out = sum(∫(1)dΓ_out) + vol_D = sum(∫(1)dΩ) + J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + C2(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out + pcfs = PDEConstrainedFunctionals(J,[C1,C2],state_map, + analytic_dJ=dJ,analytic_dC=[dC1,nothing]) + # 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; + assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg), + ls = PETScLinearSolver() + ) + # 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)) + # 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) + 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),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]) +end + +with_mpi() do distribute + mesh_partition = (5,5,5) + solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true + -ksp_converged_reason -ksp_rtol 1.0e-12" + GridapPETSc.with(args=split(solver_options)) do + main(mesh_partition,distribute) + end +end \ No newline at end of file diff --git a/scripts/_dev/Paper_scripts/jobs/elast_comp_MPI_3D.pbs b/scripts/_dev/Paper_scripts/jobs/elast_comp_MPI_3D.pbs new file mode 100644 index 00000000..38e3aaa8 --- /dev/null +++ b/scripts/_dev/Paper_scripts/jobs/elast_comp_MPI_3D.pbs @@ -0,0 +1,16 @@ +#!/bin/bash -l +#PBS -P LSTO +#PBS -N elast_comp_MPI_3D + +#PBS -l select=125:ncpus=1:mpiprocs=1:ompthreads=1:mem=8GB:cputype=6140 +#PBS -l walltime=50:00:00 +#PBS -j oe + +source $HOME/hpc-environments-main/lyra/load-ompi.sh +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/_dev/Paper_scripts/elast_comp_MPI_3D.jl diff --git a/scripts/_dev/Paper_scripts/jobs/hyperelast_comp_MPI_3D.pbs b/scripts/_dev/Paper_scripts/jobs/hyperelast_comp_MPI_3D.pbs new file mode 100644 index 00000000..a9c89845 --- /dev/null +++ b/scripts/_dev/Paper_scripts/jobs/hyperelast_comp_MPI_3D.pbs @@ -0,0 +1,16 @@ +#!/bin/bash -l +#PBS -P LSTO +#PBS -N hyperelast_comp_MPI_3D + +#PBS -l select=125:ncpus=1:mpiprocs=1:ompthreads=1:mem=8GB:cputype=6140 +#PBS -l walltime=50:00:00 +#PBS -j oe + +source $HOME/hpc-environments-main/lyra/load-ompi.sh +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/_dev/Paper_scripts/hyperelast_comp_MPI_3D.jl diff --git a/scripts/_dev/Paper_scripts/jobs/inverse_hom_MPI_3D.pbs b/scripts/_dev/Paper_scripts/jobs/inverse_hom_MPI_3D.pbs new file mode 100644 index 00000000..d1a72278 --- /dev/null +++ b/scripts/_dev/Paper_scripts/jobs/inverse_hom_MPI_3D.pbs @@ -0,0 +1,16 @@ +#!/bin/bash -l +#PBS -P LSTO +#PBS -N inverse_hom_MPI_3D + +#PBS -l select=125:ncpus=1:mpiprocs=1:ompthreads=1:mem=8GB:cputype=6140 +#PBS -l walltime=50:00:00 +#PBS -j oe + +source $HOME/hpc-environments-main/lyra/load-ompi.sh +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/_dev/Paper_scripts/inverse_hom_MPI_3D.jl diff --git a/scripts/_dev/Paper_scripts/jobs/inverter_MPI_3D.pbs b/scripts/_dev/Paper_scripts/jobs/inverter_MPI_3D.pbs new file mode 100644 index 00000000..aeeada29 --- /dev/null +++ b/scripts/_dev/Paper_scripts/jobs/inverter_MPI_3D.pbs @@ -0,0 +1,16 @@ +#!/bin/bash -l +#PBS -P LSTO +#PBS -N inverter_MPI_3D + +#PBS -l select=125:ncpus=1:mpiprocs=1:ompthreads=1:mem=8GB:cputype=6140 +#PBS -l walltime=50:00:00 +#PBS -j oe + +source $HOME/hpc-environments-main/lyra/load-ompi.sh +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/_dev/Paper_scripts/inverter_MPI_3D.jl diff --git a/scripts/_dev/Paper_scripts/jobs/therm_comp_MPI_3D.pbs b/scripts/_dev/Paper_scripts/jobs/therm_comp_MPI_3D.pbs new file mode 100644 index 00000000..795dc33d --- /dev/null +++ b/scripts/_dev/Paper_scripts/jobs/therm_comp_MPI_3D.pbs @@ -0,0 +1,17 @@ +#!/bin/bash -l + +#PBS -P LSTO +#PBS -N therm_comp_MPI_3D + +#PBS -l select=125:ncpus=1:mpiprocs=1:ompthreads=1:mem=8GB:cputype=6140 +#PBS -l walltime=23:00:00 +#PBS -j oe + +source $HOME/hpc-environments-main/lyra/load-ompi.sh +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/_dev/Paper_scripts/therm_comp_MPI_3D.jl \ No newline at end of file diff --git a/scripts/_dev/Added_for_paper/thermal_MPI_PETSc.jl b/scripts/_dev/Paper_scripts/therm_comp_MPI.jl similarity index 86% rename from scripts/_dev/Added_for_paper/thermal_MPI_PETSc.jl rename to scripts/_dev/Paper_scripts/therm_comp_MPI.jl index ef08aa72..a0b1afeb 100644 --- a/scripts/_dev/Added_for_paper/thermal_MPI_PETSc.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_MPI.jl @@ -5,7 +5,7 @@ function main(mesh_partition,distribute) ranks = distribute(LinearIndices((prod(mesh_partition),))) # FE parameters order = 1 # Finite element order - xmax = ymax = 1.0 # Domain size + xmax,ymax = (1.0,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 @@ -25,7 +25,7 @@ function main(mesh_partition,distribute) 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_MPI_PETSc/" # Output path + path = "./results/therm_comp_MPI/" # Output path i_am_main(ranks) && mkpath(path) # Create path # Model model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); @@ -64,10 +64,10 @@ function main(mesh_partition,distribute) 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]) + C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + pcfs = PDEConstrainedFunctionals(J,[C1],state_map, + analytic_dJ=dJ,analytic_dC=[dC1]) # Velocity extension α = 4*maximum(get_el_Δ(model)) a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ @@ -80,12 +80,11 @@ function main(mesh_partition,distribute) 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]) + optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=i_am_main(ranks)) # 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()) + iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data) write_history(path*"/history.txt",get_history(optimiser);ranks) end # Final structure diff --git a/scripts/_dev/Added_for_paper/thermal_MPI_PETSc_3D.jl b/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl similarity index 87% rename from scripts/_dev/Added_for_paper/thermal_MPI_PETSc_3D.jl rename to scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl index 90addf1d..0fb2889c 100644 --- a/scripts/_dev/Added_for_paper/thermal_MPI_PETSc_3D.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_MPI_3D.jl @@ -5,8 +5,8 @@ function main(mesh_partition,distribute) 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 + xmax,ymax,zmax = (1.0,1.0,1.0) # Domain size + dom = (0,xmax,0,ymax,0,zmax) # Bounding domain el_size = (150,150,150) # Mesh partition size prop_Γ_N = 0.2 # Γ_N size parameter prop_Γ_D = 0.2 # Γ_D size parameter @@ -27,7 +27,7 @@ function main(mesh_partition,distribute) 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_MPI_PETSc_3D/" # Output path + path = "./results/therm_comp_MPI_3D/" # Output path i_am_main(ranks) && mkpath(path) # Create path # Model model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); @@ -66,10 +66,10 @@ function main(mesh_partition,distribute) 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]) + C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + pcfs = PDEConstrainedFunctionals(J,[C1],state_map, + analytic_dJ=dJ,analytic_dC=[dC1]) # Velocity extension α = 4*maximum(get_el_Δ(model)) a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ @@ -82,12 +82,12 @@ function main(mesh_partition,distribute) 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]) + optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + verbose=i_am_main(ranks)) # 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()) + iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data) write_history(path*"/history.txt",get_history(optimiser);ranks) end # Final structure diff --git a/scripts/_dev/Added_for_paper/thermal_serial.jl b/scripts/_dev/Paper_scripts/therm_comp_serial.jl similarity index 93% rename from scripts/_dev/Added_for_paper/thermal_serial.jl rename to scripts/_dev/Paper_scripts/therm_comp_serial.jl index 115e1b5e..055f5ad9 100644 --- a/scripts/_dev/Added_for_paper/thermal_serial.jl +++ b/scripts/_dev/Paper_scripts/therm_comp_serial.jl @@ -23,7 +23,7 @@ function main() 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/" # Output path + path = "./results/therm_comp_serial/" # Output path mkpath(path) # Create path # Model model = CartesianDiscreteModel(dom,el_size); @@ -65,8 +65,8 @@ function main() 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]) + optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit, + verbose=true) # Solve for (it,uh,φh) in optimiser data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]