diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl index c0c0dbea..c98cca21 100644 --- a/src/GridapTopOpt.jl +++ b/src/GridapTopOpt.jl @@ -61,8 +61,6 @@ export reinit! include("Solvers.jl") export ElasticitySolver -export BlockDiagonalPreconditioner -export MUMPSSolver include("VelocityExtension.jl") export VelocityExtension diff --git a/src/Solvers.jl b/src/Solvers.jl index a81ab340..b1790443 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -1,4 +1,4 @@ -## TODO: @Jordi, are these going into GridapSolvers or should I write documentation? +## This will be moved into GridapSolvers in a future release struct ElasticitySolver{A} <: LinearSolver space ::A @@ -31,7 +31,7 @@ function get_dof_coordinates(space::GridapDistributed.DistributedSingleFieldFESp owner = part_id(dof_indices) own_indices = OwnIndices(ngdofs,owner,own_to_global(dof_indices)) ghost_indices = GhostIndices(ngdofs,Int64[],Int32[]) # We only consider owned dofs - OwnAndGhostIndices(own_indices,ghost_indices) + OwnAndGhostIndices(own_indices,ghost_indices) end return PVector(coords,indices) end @@ -126,7 +126,7 @@ _num_dims(space::GridapDistributed.DistributedSingleFieldFESpace) = getany(map(_ function Gridap.Algebra.numerical_setup(ss::ElasticitySymbolicSetup,_A::PSparseMatrix) s = ss.solver - # Create ns + # Create ns A = convert(PETScMatrix,_A) X = convert(PETScVector,allocate_in_domain(_A)) B = convert(PETScVector,allocate_in_domain(_A)) @@ -139,7 +139,7 @@ function Gridap.Algebra.numerical_setup(ss::ElasticitySymbolicSetup,_A::PSparseM # Create matrix nullspace @check_error_code GridapPETSc.PETSC.MatNullSpaceCreateRigidBody(dof_coords.vec[],ns.null) @check_error_code GridapPETSc.PETSC.MatSetNearNullSpace(ns.A.mat[],ns.null[]) - + # Setup solver and preconditioner @check_error_code GridapPETSc.PETSC.KSPCreate(ns.A.comm,ns.ksp) @check_error_code GridapPETSc.PETSC.KSPSetOperators(ns.ksp[],ns.A.mat[],ns.A.mat[]) @@ -162,81 +162,4 @@ function Algebra.solve!(x::AbstractVector{PetscScalar},ns::ElasticityNumericalSe @check_error_code GridapPETSc.PETSC.KSPSolve(ns.ksp[],B.vec[],X.vec[]) copy!(x,X) return x -end - -# MUMPS solver - -function mumps_ksp_setup(ksp) - pc = Ref{GridapPETSc.PETSC.PC}() - mumpsmat = Ref{GridapPETSc.PETSC.Mat}() - @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL) - @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPPREONLY) - @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc) - @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCLU) - @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[],GridapPETSc.PETSC.MATSOLVERMUMPS) - @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[]) - @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat) - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 4, 1) # Level of printing - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 7, 0) # Perm type -# @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 14, 1000) - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2) # Seq or parallel analysis - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 2) # Parallel ordering - @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6) # Absolute pivoting threshold -end - -function MUMPSSolver() - return PETScLinearSolver(mumps_ksp_setup) -end - -# Block diagonal preconditioner - -struct BlockDiagonalPreconditioner{N,A} <: Gridap.Algebra.LinearSolver - solvers :: A - function BlockDiagonalPreconditioner(solvers::AbstractArray{<:Gridap.Algebra.LinearSolver}) - N = length(solvers) - A = typeof(solvers) - return new{N,A}(solvers) - end -end - -struct BlockDiagonalPreconditionerSS{A,B} <: Gridap.Algebra.SymbolicSetup - solver :: A - block_ss :: B -end - -function Gridap.Algebra.symbolic_setup(solver::BlockDiagonalPreconditioner,mat::AbstractBlockMatrix) - mat_blocks = diag(blocks(mat)) - block_ss = map(symbolic_setup,solver.solvers,mat_blocks) - return BlockDiagonalPreconditionerSS(solver,block_ss) -end - -struct BlockDiagonalPreconditionerNS{A,B} <: Gridap.Algebra.NumericalSetup - solver :: A - block_ns :: B -end - -function Gridap.Algebra.numerical_setup(ss::BlockDiagonalPreconditionerSS,mat::AbstractBlockMatrix) - solver = ss.solver - mat_blocks = diag(blocks(mat)) - block_ns = map(numerical_setup,ss.block_ss,mat_blocks) - return BlockDiagonalPreconditionerNS(solver,block_ns) -end - -function Gridap.Algebra.numerical_setup!(ns::BlockDiagonalPreconditionerNS,mat::AbstractBlockMatrix) - mat_blocks = diag(blocks(mat)) - map(numerical_setup!,ns.block_ns,mat_blocks) -end - -function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::BlockDiagonalPreconditionerNS,b::AbstractBlockVector) - @check blocklength(x) == blocklength(b) == length(ns.block_ns) - for (iB,bns) in enumerate(ns.block_ns) - xi = x[Block(iB)] - bi = b[Block(iB)] - solve!(xi,bns,bi) - end - return x -end - -function LinearAlgebra.ldiv!(x,ns::BlockDiagonalPreconditionerNS,b) - solve!(x,ns,b) -end +end \ No newline at end of file diff --git a/test/mpi/InverseHomogenisationALMTests.jl b/test/mpi/InverseHomogenisationALMTests.jl new file mode 100644 index 00000000..dcac0648 --- /dev/null +++ b/test/mpi/InverseHomogenisationALMTests.jl @@ -0,0 +1,108 @@ +module InverseHomogenisationALMTestsMPI +using Test + +using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers, + PartitionedArrays, GridapTopOpt, SparseMatricesCSR + +""" + (MPI) Maximum bulk modulus inverse homogenisation with augmented Lagrangian method in 2D. + + Optimisation problem: + Min J(Ω) = -κ(Ω) + Ω + s.t., Vol(Ω) = vf, + ⎡For unique εᴹᵢ, find uᵢ∈V=H¹ₚₑᵣ(Ω)ᵈ, + ⎣∫ ∑ᵢ C ⊙ ε(uᵢ) ⊙ ε(vᵢ) dΩ = ∫ -∑ᵢ C ⊙ ε⁰ᵢ ⊙ ε(vᵢ) dΩ, ∀v∈V. +""" +function main(distribute,mesh_partition;AD) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + order = 1 + xmax,ymax=(1.0,1.0) + dom = (0,xmax,0,ymax) + el_size = (10,10) + γ = 0.05 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + C = isotropic_elast_tensor(2,1.,0.3) + η_coeff = 2 + α_coeff = 4max_steps*γ + vf = 0.5 + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size,isperiodic=(true,true)) + el_Δ = get_el_Δ(model) + f_Γ_D(x) = iszero(x) + update_labels!(1,model,f_Γ_D,"origin") + + ## Triangulations and measures + Ω = Triangulation(model) + dΩ = Measure(Ω,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=["origin"]) + U = TrialFESpace(V,VectorValue(0.0,0.0)) + V_reg = V_φ = TestFESpace(model,reffe_scalar) + U_reg = TrialFESpace(V_reg) + + ## Create FE functions + lsf_fn = x->max(initial_lsf(2,0.4)(x),initial_lsf(2,0.4;b=VectorValue(0,0.5))(x)) + φh = interpolate(lsf_fn,V_φ) + + ## Interpolation and weak form + interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + + εᴹ = (TensorValue(1.,0.,0.,0.), # ϵᵢⱼ⁽¹¹⁾≡ϵᵢⱼ⁽¹⁾ + TensorValue(0.,0.,0.,1.), # ϵᵢⱼ⁽²²⁾≡ϵᵢⱼ⁽²⁾ + TensorValue(0.,1/2,1/2,0.)) # ϵᵢⱼ⁽¹²⁾≡ϵᵢⱼ⁽³⁾ + + a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ + l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3] + + ## 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Ω + κ(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Ω + + ## Finite difference solver and level set function + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + + ## Setup solver and FE operators + state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ) + pcfs = if AD + PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) + else + PDEConstrainedFunctionals(κ,[Vol],state_map) + end + + ## 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,verbose=i_am_main(ranks),constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +with_mpi() do distribute + @test main(distribute,(2,2);AD=true) + @test main(distribute,(2,2);AD=false) +end + + +end # module \ No newline at end of file diff --git a/test/mpi/InverterHPMTests.jl b/test/mpi/InverterHPMTests.jl new file mode 100644 index 00000000..055a9f20 --- /dev/null +++ b/test/mpi/InverterHPMTests.jl @@ -0,0 +1,118 @@ +module InverterHPMTestsMPI +using Test + +using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers, + PartitionedArrays, GridapTopOpt, SparseMatricesCSR + +""" + (MPI) Inverter mechanism with Hilbertian projection method in 2D. + + Optimisation problem: + Min J(Ω) = ηᵢₙ*∫ u⋅e₁ dΓᵢₙ/Vol(Γᵢₙ) + Ω + s.t., Vol(Ω) = vf, + C(Ω) = 0, + ⎡u∈V=H¹(Ω;u(Γ_D)=0)ᵈ, + ⎣∫ C ⊙ ε(u) ⊙ ε(v) dΩ + ∫ kₛv⋅u dΓₒᵤₜ = ∫ v⋅g dΓᵢₙ , ∀v∈V. + + where C(Ω) = ∫ -u⋅e₁-δₓ dΓₒᵤₜ/Vol(Γₒᵤₜ). We assume symmetry in the problem to aid + convergence. +""" +function main(distribute,mesh_partition) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + order = 1 + dom = (0,1,0,0.5) + el_size = (10,10) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + C = isotropic_elast_tensor(2,1.0,0.3) + η_coeff = 2 + α_coeff = 4max_steps*γ + vf = 0.4 + δₓ = 0.2 + ks = 0.1 + g = VectorValue(0.5,0) + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size) + el_Δ = get_el_Δ(model) + f_Γ_in(x) = (x[1] ≈ 0.0) && (x[2] <= 0.03 + eps()) + f_Γ_out(x) = (x[1] ≈ 1.0) && (x[2] <= 0.07 + eps()) + f_Γ_D(x) = (x[1] ≈ 0.0) && (x[2] >= 0.4) + f_Γ_D2(x) = (x[2] ≈ 0.0) + update_labels!(1,model,f_Γ_in,"Gamma_in") + update_labels!(2,model,f_Γ_out,"Gamma_out") + update_labels!(3,model,f_Γ_D,"Gamma_D") + update_labels!(4,model,f_Γ_D2,"SymLine") + + ## Triangulations and measures + Ω = Triangulation(model) + Γ_in = BoundaryTriangulation(model,tags="Gamma_in") + Γ_out = BoundaryTriangulation(model,tags="Gamma_out") + dΩ = Measure(Ω,2order) + dΓ_in = Measure(Γ_in,2order) + dΓ_out = Measure(Γ_out,2order) + vol_D = sum(∫(1)dΩ) + vol_Γ_in = sum(∫(1)dΓ_in) + vol_Γ_out = sum(∫(1)dΓ_out) + + ## Spaces + reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D","SymLine"], + dirichlet_masks=[(true,true),(false,true)]) + U = TrialFESpace(V,[VectorValue(0.0,0.0),VectorValue(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,0]) + + ## Create FE functions + lsf_fn(x) = min(max(initial_lsf(6,0.2)(x),-sqrt((x[1]-1)^2+(x[2]-0.5)^2)+0.2), + sqrt((x[1])^2+(x[2]-0.5)^2)-0.1) + φh = interpolate(lsf_fn,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Γ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in + + ## Optimisation functionals + e₁ = VectorValue(1,0) + J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out + + ## Finite difference solver and level set function + ls_evo = HamiltonJacobiEvolution(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Γ_in,dΓ_out) + pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing]) + + ## 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 = HilbertianProjection(pcfs,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),debug=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +with_mpi() do distribute + @test main(distribute,(2,2)) +end + +end # module \ No newline at end of file diff --git a/test/mpi/NeohookAnalyticJacALMTests.jl b/test/mpi/NeohookAnalyticJacALMTests.jl new file mode 100644 index 00000000..df8041ac --- /dev/null +++ b/test/mpi/NeohookAnalyticJacALMTests.jl @@ -0,0 +1,125 @@ +module NeohookAnalyticJacALMTestsMPI +using Test + +using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers, + PartitionedArrays, GridapTopOpt, SparseMatricesCSR + +""" + (MPI) Minimum hyperelastic (neohookean) compliance with augmented Lagrangian method in 2D. +""" +function main(distribute,mesh_partition) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + order = 1 + xmax=ymax=1.0 + prop_Γ_N = 0.2 + dom = (0,xmax,0,ymax) + el_size = (10,10) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + η_coeff = 2 + α_coeff = 4max_steps*γ + vf = 0.6 + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,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 - ν)) + g = VectorValue(0,-20) + + ## 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 + + function dS(∇du,∇u) + Cinv = inv(C(F(∇u))) + _dE = dE(∇du,∇u) + λ*(Cinv⊙_dE)*Cinv + 2*(μ-λ*log(J(F(∇u))))*Cinv⋅_dE⋅(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 + jac_mat(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ + jac_geo(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ + jac(u,du,v,φ,dΩ,dΓ_N) = jac_mat(u,du,v,φ,dΩ,dΓ_N) + jac_geo(u,du,v,φ,dΩ,dΓ_N) + + ## Optimisation functionals + J(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 + ls_evo = HamiltonJacobiEvolution(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;jac) + pcfs = PDEConstrainedFunctionals(J,[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,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +with_mpi() do distribute + @test main(distribute,(2,2)) +end + +end # module \ No newline at end of file diff --git a/test/mpi/NonlinearThermalComplianceALMTests.jl b/test/mpi/NonlinearThermalComplianceALMTests.jl new file mode 100644 index 00000000..1b63a033 --- /dev/null +++ b/test/mpi/NonlinearThermalComplianceALMTests.jl @@ -0,0 +1,105 @@ +module NonlinearThermalComplianceALMTestsMPI +using Test + +using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers, + PartitionedArrays, GridapTopOpt, SparseMatricesCSR + +""" + (MPI) Minimum thermal compliance with augmented Lagrangian method in 2D with nonlinear diffusivity. + + Optimisation problem: + Min J(Ω) = ∫ κ(u)*∇(u)⋅∇(u) dΩ + Ω + s.t., Vol(Ω) = vf, + ⎡u∈V=H¹(Ω;u(Γ_D)=0), + ⎣∫ κ(u)*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. + + In this example κ(u) = κ0*(exp(ξ*u)) +""" +function main(distribute,mesh_partition) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + order = 1 + xmax=ymax=1.0 + prop_Γ_N = 0.2 + prop_Γ_D = 0.2 + dom = (0,xmax,0,ymax) + el_size = (10,10) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + η_coeff = 2 + α_coeff = 4max_steps*γ + vf = 0.4 + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); + el_Δ = get_el_Δ(model) + f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= ymax*prop_Γ_D + eps() || + x[2] >= ymax-ymax*prop_Γ_D - eps())); + 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_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) + + ## 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.ρ + + κ0 = 1 + ξ = -1 + κ(u) = κ0*(exp(ξ*u)) + res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N + + ## Optimisation functionals + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(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 + ls_evo = HamiltonJacobiEvolution(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(J,[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,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +with_mpi() do distribute + @test main(distribute,(2,2)) +end + +end # module \ No newline at end of file diff --git a/test/mpi/PZMultiFieldRepeatingStateTests.jl b/test/mpi/PZMultiFieldRepeatingStateTests.jl new file mode 100644 index 00000000..71c2e623 --- /dev/null +++ b/test/mpi/PZMultiFieldRepeatingStateTests.jl @@ -0,0 +1,202 @@ +module PZMultiFieldRepeatingStateTestsMPI +using Test + +using Gridap, GridapDistributed, GridapPETSc, GridapSolvers, + PartitionedArrays, GridapTopOpt, SparseMatricesCSR + +using Gridap.TensorValues, Gridap.Helpers, Gridap.MultiField + +## Parameters +function main(distribute,mesh_partition;AD,use_mfs) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + el_size = (10,10) + order = 1 + xmax,ymax=(1.0,1.0) + dom = (0,xmax,0,ymax) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + η_coeff = 2 + α_coeff = 4*max_steps*γ + vf = 0.5 + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size,isperiodic=(true,true)) + el_Δ = get_el_Δ(model) + f_Γ_D(x) = iszero(x) + update_labels!(1,model,f_Γ_D,"origin") + + ## Triangulations and measures + Ω = Triangulation(model) + dΩ = Measure(Ω,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;conformity=:H1,dirichlet_tags=["origin"]) + U = TrialFESpace(V,VectorValue(0.0,0.0)) + Q = TestFESpace(model,reffe_scalar;conformity=:H1,dirichlet_tags=["origin"]) + P = TrialFESpace(Q,0) + if use_mfs + mfs = BlockMultiFieldStyle() + UP = MultiFieldFESpace([U,P];style=mfs) + VQ = MultiFieldFESpace([V,Q];style=mfs) + else + UP = MultiFieldFESpace([U,P]) + VQ = MultiFieldFESpace([V,Q]) + end + + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar) + U_reg = TrialFESpace(V_reg) + + ## Create FE functions + φh = interpolate(initial_lsf(2,0.2),V_φ) + + ## Interpolation and weak form + interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ),ϵ=10^-9) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + + ## Material tensors + C, e, κ = PZT5A_2D(); + k0 = norm(C.data,Inf); + α0 = norm(e.data,Inf); + β0 = norm(κ.data,Inf); + γ0 = β0*k0/α0^2; + + ## Weak forms + εᴹ = (SymTensorValue(1.0,0.0,0.0), + SymTensorValue(0.0,0.0,1.0), + SymTensorValue(0.0,1/2,0)) + Eⁱ = (VectorValue(1.0,0.0,), + VectorValue(0.0,1.0)) + + a((u,ϕ),(v,q),φ,dΩ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) - + 1/α0*((-∇(ϕ) ⋅ e) ⊙ ε(v)) + + -1/α0*((e ⋅² ε(u)) ⋅ -∇(q)) + + -γ0/β0*((κ ⋅ -∇(ϕ)) ⋅ -∇(q))) )dΩ; + + l_ε = [((v,q),φ,dΩ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3]; + l_E = [((v,q),φ,dΩ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2]; + l = [l_ε; l_E] + + function Cᴴ(r,s,uϕ,φ,dΩ) + u_s = uϕ[2s-1]; ϕ_s = uϕ[2s] + ∫(1/k0 * (I ∘ φ) * (((C ⊙ (1/k0*ε(u_s) + εᴹ[s])) ⊙ εᴹ[r]) - ((-1/α0*∇(ϕ_s) ⋅ e) ⊙ εᴹ[r])))dΩ; + end + + function DCᴴ(r,s,q,uϕ,φ,dΩ) + u_r = uϕ[2r-1]; ϕ_r = uϕ[2r] + u_s = uϕ[2s-1]; ϕ_s = uϕ[2s] + ∫(- 1/k0 * q * ( + (C ⊙ (1/k0*ε(u_s) + εᴹ[s])) ⊙ (1/k0*ε(u_r) + εᴹ[r]) - + (-1/α0*∇(ϕ_s) ⋅ e) ⊙ (1/k0*ε(u_r) + εᴹ[r]) - + (e ⋅² (1/k0*ε(u_s) + εᴹ[s])) ⋅ (-1/α0*∇(ϕ_r)) - + (κ ⋅ (-1/α0*∇(ϕ_s))) ⋅ (-1/α0*∇(ϕ_r)) + ) * (DH ∘ φ) * (norm ∘ ∇(φ)) + )dΩ; + end + + Bᴴ(uϕ,φ,dΩ) = 1/4*(Cᴴ(1,1,uϕ,φ,dΩ)+Cᴴ(2,2,uϕ,φ,dΩ)+2*Cᴴ(1,2,uϕ,φ,dΩ)) + DBᴴ(q,uϕ,φ,dΩ) = 1/4*(DCᴴ(1,1,q,uϕ,φ,dΩ)+DCᴴ(2,2,q,uϕ,φ,dΩ)+2*DCᴴ(1,2,q,uϕ,φ,dΩ)) + + J(uϕ,φ,dΩ) = -1*Bᴴ(uϕ,φ,dΩ) + DJ(q,uϕ,φ,dΩ) = -1*DBᴴ(q,uϕ,φ,dΩ) + C1(uϕ,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + DC1(q,uϕ,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + + ## Finite difference solver and level set function + stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + + ## Setup solver and FE operators + state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh,dΩ) + pcfs = if AD + PDEConstrainedFunctionals(J,[C1],state_map;analytic_dJ=DJ,analytic_dC=[DC1]) + else + PDEConstrainedFunctionals(J,[C1],state_map) + end + + ## 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 = HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=i_am_main(ranks)) + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +function PZT5A_2D() + ε_0 = 8.854e-12; + C_Voigt = [12.0400e10 7.51000e10 0.0 + 7.51000e10 11.0900e10 0.0 + 0.0 0.0 2.1000e10] + e_Voigt = [0.0 0.0 12.30000 + -5.40000 15.80000 0.0] + K_Voigt = [540*ε_0 0 + 0 830*ε_0] + C = voigt2tensor4(C_Voigt) + e = voigt2tensor3(e_Voigt) + κ = voigt2tensor2(K_Voigt) + C,e,κ +end + +""" + Given a material constant given in Voigt notation, + return a SymFourthOrderTensorValue using ordering from Gridap +""" +function voigt2tensor4(A::Array{M,2}) where M + if isequal(size(A),(3,3)) + return SymFourthOrderTensorValue(A[1,1], A[3,1], A[2,1], + A[1,3], A[3,3], A[2,3], + A[1,2], A[3,2], A[2,2]) + elseif isequal(size(A),(6,6)) + return SymFourthOrderTensorValue(A[1,1], A[6,1], A[5,1], A[2,1], A[4,1], A[3,1], + A[1,6], A[6,6], A[5,6], A[2,6], A[4,6], A[3,6], + A[1,5], A[6,5], A[5,5], A[2,5], A[4,5], A[3,5], + A[1,2], A[6,2], A[5,2], A[2,2], A[4,2], A[3,2], + A[1,4], A[6,4], A[5,4], A[2,4], A[4,4], A[3,4], + A[1,3], A[6,3], A[5,3], A[2,3], A[4,3], A[3,3]) + else + @notimplemented + end +end + +""" + Given a material constant given in Voigt notation, + return a ThirdOrderTensorValue using ordering from Gridap +""" +function voigt2tensor3(A::Array{M,2}) where M + if isequal(size(A),(2,3)) + return ThirdOrderTensorValue(A[1,1], A[2,1], A[1,3], A[2,3], A[1,3], A[2,3], A[1,2], A[2,2]) + elseif isequal(size(A),(3,6)) + return ThirdOrderTensorValue( + A[1,1], A[2,1], A[3,1], A[1,6], A[2,6], A[3,6], A[1,5], A[2,5], A[3,5], + A[1,6], A[2,6], A[3,6], A[1,2], A[2,2], A[3,2], A[1,4], A[2,4], A[3,4], + A[1,5], A[2,5], A[3,5], A[1,4], A[2,4], A[3,4], A[1,3], A[2,3], A[3,3]) + else + @notimplemented + end +end + +""" + Given a material constant given in Voigt notation, + return a SymTensorValue using ordering from Gridap +""" +function voigt2tensor2(A::Array{M,2}) where M + return TensorValue(A) +end + +# Test that these run successfully +with_mpi() do distribute + @test main(distribute,(2,2);AD=true) + @test main(distribute,(2,2);AD=false) + @test main(distribute,(2,2);AD=true,use_mfs=true) + @test main(distribute,(2,2);AD=false,use_mfs=true) +end + +end # module \ No newline at end of file diff --git a/test/mpi/ThermalComplianceALMTests.jl b/test/mpi/ThermalComplianceALMTests.jl new file mode 100644 index 00000000..dd4863d8 --- /dev/null +++ b/test/mpi/ThermalComplianceALMTests.jl @@ -0,0 +1,199 @@ +module ThermalComplianceALMTestsMPI +using Test + +using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers, + PartitionedArrays, GridapTopOpt, SparseMatricesCSR + +""" + (MPI) Minimum thermal compliance with augmented Lagrangian method in 2D. + + Optimisation problem: + Min J(Ω) = ∫ κ*∇(u)⋅∇(u) dΩ + Ω + s.t., Vol(Ω) = vf, + ⎡u∈V=H¹(Ω;u(Γ_D)=0), + ⎣∫ κ*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. +""" +function main(distribute,mesh_partition;order,AD) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + xmax = ymax = 1.0 + prop_Γ_N = 0.2 + prop_Γ_D = 0.2 + dom = (0,xmax,0,ymax) + el_size = (10,10) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5*order^2)/minimum(el_size) + κ = 1 + vf = 0.4 + η_coeff = 2 + α_coeff = 4max_steps*γ + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); + el_Δ = get_el_Δ(model) + f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= ymax*prop_Γ_D + eps() || + x[2] >= ymax-ymax*prop_Γ_D - eps())) + 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_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) + + ## 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 ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N + + ## Optimisation functionals + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(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 + ls_evo = HamiltonJacobiEvolution(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 = if AD + PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) + else + PDEConstrainedFunctionals(J,[Vol],state_map) + end + + ## 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,verbose=i_am_main(ranks),constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +""" + (MPI) Minimum thermal compliance with augmented Lagrangian method in 3D. + + Optimisation problem: + Min J(Ω) = ∫ κ*∇(u)⋅∇(u) dΩ + Ω + s.t., Vol(Ω) = vf, + ⎡u∈V=H¹(Ω;u(Γ_D)=0), + ⎣∫ κ*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. +""" +function main_3d(ranks,mesh_partition,;order) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + xmax = ymax = zmax = 1.0 + prop_Γ_N = 0.2 + prop_Γ_D = 0.2 + dom = (0,xmax,0,ymax,0,zmax) + el_size = (20,20,20) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5*order^2)/minimum(el_size) + κ = 1 + vf = 0.4 + η_coeff = 2 + α_coeff = 4max_steps*γ + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); + el_Δ = get_el_Δ(model) + f_Γ_D(x) = (x[1] ≈ 0.0) && (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()) + f_Γ_N(x) = (x[1] ≈ xmax) && (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()) + 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) + + ## 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 ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N + + ## Optimisation functionals + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(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 + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,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,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +with_mpi() do distribute + @test main(distribute,(2,2);order=1,AD=true) + @test main(distribute,(2,2);order=2,AD=true) + @test main(distribute,(2,2);order=1,AD=false) + @test main_3d(distribute,(2,2);order=1) +end + +end # module \ No newline at end of file diff --git a/test/mpi/ThermalComplianceALMTestsPETSc.jl b/test/mpi/ThermalComplianceALMTestsPETSc.jl new file mode 100644 index 00000000..f04091bf --- /dev/null +++ b/test/mpi/ThermalComplianceALMTestsPETSc.jl @@ -0,0 +1,229 @@ +module ThermalComplianceALMTestsMPIPETSc +using Test + +using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers, + PartitionedArrays, GridapTopOpt, SparseMatricesCSR + +""" + (MPI) Minimum thermal compliance with augmented Lagrangian method in 2D. + + Optimisation problem: + Min J(Ω) = ∫ κ*∇(u)⋅∇(u) dΩ + Ω + s.t., Vol(Ω) = vf, + ⎡u∈V=H¹(Ω;u(Γ_D)=0), + ⎣∫ κ*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. +""" +function main(distribute,mesh_partition;order,AD) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + xmax = ymax = 1.0 + prop_Γ_N = 0.2 + prop_Γ_D = 0.2 + dom = (0,xmax,0,ymax) + el_size = (10,10) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5*order^2)/minimum(el_size) + κ = 1 + vf = 0.4 + η_coeff = 2 + α_coeff = 4max_steps*γ + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); + el_Δ = get_el_Δ(model) + f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= ymax*prop_Γ_D + eps() || + x[2] >= ymax-ymax*prop_Γ_D - eps())) + 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_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) + + ## 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 ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N + + ## Optimisation functionals + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(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 + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + + ## Setup solver and FE operators + 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 + ) + pcfs = if AD + PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) + else + PDEConstrainedFunctionals(J,[Vol],state_map) + end + + ## 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; + assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg), + ls = solver + ) + + ## Optimiser + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +""" + (MPI) Minimum thermal compliance with augmented Lagrangian method in 3D. + + Optimisation problem: + Min J(Ω) = ∫ κ*∇(u)⋅∇(u) dΩ + Ω + s.t., Vol(Ω) = vf, + ⎡u∈V=H¹(Ω;u(Γ_D)=0), + ⎣∫ κ*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. +""" +function main_3d(ranks,mesh_partition,;order) + ranks = distribute(LinearIndices((prod(mesh_partition),))) + ## Parameters + xmax = ymax = zmax = 1.0 + prop_Γ_N = 0.2 + prop_Γ_D = 0.2 + dom = (0,xmax,0,ymax,0,zmax) + el_size = (20,20,20) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5*order^2)/minimum(el_size) + κ = 1 + vf = 0.4 + η_coeff = 2 + α_coeff = 4max_steps*γ + + ## FE Setup + model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size); + el_Δ = get_el_Δ(model) + f_Γ_D(x) = (x[1] ≈ 0.0) && (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()) + f_Γ_N(x) = (x[1] ≈ xmax) && (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()) + 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) + + ## 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 ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N + + ## Optimisation functionals + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(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 + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + + ## Setup solver and FE operators + 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 + ) + 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; + assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg), + ls = solver + ) + + ## Optimiser + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +with_mpi() do distribute + 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 + @test main(distribute,(2,2);order=1,AD=false) + @test main_3d(distribute,(2,2);order=1) + end +end + +end # module \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 64966857..cf6e0ec3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,2 +1,3 @@ using Test -include("seq/runtests.jl") \ No newline at end of file +include("seq/runtests.jl") +include("mpi/runtests.jl") \ No newline at end of file