-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Added MPI tests and cleaned up Solvers.jl
- Loading branch information
Showing
10 changed files
with
1,093 additions
and
85 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
Oops, something went wrong.