Merge main into src_refactor_and_docs

.gitignore

@@ -19,8 +19,6 @@ deps/src/
 # Build artifacts for creating documentation generated by the Documenter package
 docs/build/
 docs/site/
-Results/*
-results/*
 
 # File generated by Pkg, the package manager, based on a corresponding Project.toml
 # It records a fixed state of all packages used by the project. As such, it should not be

results/.gitignore

@@ -0,0 +1,2 @@
+*
+!.gitignore

scripts/Benchmarks/benchmark.jl

@@ -439,7 +439,7 @@ with_mpi() do distribute
       bopt0 = i_am_main(ranks) && zeros(NREPS)
     end
     if occursin("bopt1",BMARK_TYPE)
-      bopt1 = benchmark_optimizer(optim, 1, ranks; nreps=NREPS)
+      bopt1 = benchmark_single_iteration(optim, ranks; nreps = NREPS)
     else
       bopt1 = i_am_main(ranks) && zeros(NREPS)
     end

scripts/MPI/3d_elastic_compliance_ALM.jl

@@ -64,7 +64,6 @@ function main(mesh_partition,distribute,el_size)
 
   a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
   l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N
-  res(u,v,φ,dΩ,dΓ_N) = a(u,v,φ,dΩ,dΓ_N) - l(v,φ,dΩ,dΓ_N)
 
   ## Optimisation functionals
   J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ

scripts/MPI/3d_inverse_homenisation_ALM.jl

@@ -44,7 +44,7 @@ function main(mesh_partition,distribute,el_size)
   reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
   reffe_scalar = ReferenceFE(lagrangian,Float64,order)
   V = TestFESpace(model,reffe;dirichlet_tags=["origin"])
-  U = TrialFESpace(_V,VectorValue(0.0,0.0,0.0))
+  U = TrialFESpace(V,VectorValue(0.0,0.0,0.0))
   V_reg = V_φ = TestFESpace(model,reffe_scalar)
   U_reg = TrialFESpace(V_reg)
 
@@ -67,16 +67,12 @@ function main(mesh_partition,distribute,el_size)
   l = [(v,φ,dΩ) -> ∫(-(I ∘ φ) * C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:6]
 
   ## Optimisation functionals
-  _C(C,ε_p,ε_q) = C ⊙ ε_p ⊙ ε_q;
-
-  _K(C,u,εᴹ) = (_C(C,ε(u[1])+εᴹ[1],εᴹ[1]) + _C(C,ε(u[2])+εᴹ[2],εᴹ[2]) + _C(C,ε(u[3])+εᴹ[3],εᴹ[3]) + 
-              2(_C(C,ε(u[1])+εᴹ[1],εᴹ[2]) + _C(C,ε(u[1])+εᴹ[1],εᴹ[3]) + _C(C,ε(u[2])+εᴹ[2],εᴹ[3])))/9 
-
-  _v_K(C,u,εᴹ) = (_C(C,ε(u[1])+εᴹ[1],ε(u[1])+εᴹ[1]) + _C(C,ε(u[2])+εᴹ[2],ε(u[2])+εᴹ[2]) + _C(C,ε(u[3])+εᴹ[3],ε(u[3])+εᴹ[3]) + 
-                2(_C(C,ε(u[1])+εᴹ[1],ε(u[2])+εᴹ[2]) + _C(C,ε(u[1])+εᴹ[1],ε(u[3])+εᴹ[3]) + _C(C,ε(u[2])+εᴹ[2],ε(u[3])+εᴹ[3])))/9 
-
-  J(u,φ,dΩ) = ∫(-(I ∘ φ)*_K(C,u,εᴹ))dΩ
-  dJ(q,u,φ,dΩ) = ∫(_v_K(C,u,εᴹ)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ,dΩ) = -1/9*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+Cᴴ(3,3,u,φ,dΩ)+
+    2*(Cᴴ(1,2,u,φ,dΩ)+Cᴴ(1,3,u,φ,dΩ)+Cᴴ(2,3,u,φ,dΩ)))
+  dκ(q,u,φ,dΩ) = -1/9*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+dCᴴ(3,3,q,u,φ,dΩ)+
+    2*(dCᴴ(1,2,q,u,φ,dΩ)+dCᴴ(1,3,q,u,φ,dΩ)+dCᴴ(2,3,q,u,φ,dΩ)))
   Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
   dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
@@ -95,7 +91,7 @@ function main(mesh_partition,distribute,el_size)
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
     ls = solver, adjoint_ls = solver
   )
-  pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dJ=dJ,analytic_dC=[dVol])
+  pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
   α = α_coeff*maximum(el_Δ)

scripts/MPI/3d_inverter_ALM.jl

@@ -117,9 +117,8 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  ls_enabled=false # Setting to true will use a line search instead of oscillation detection
-  optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,α_min=0.4,ls_enabled,
-    verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out])
+  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+    γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out])
   for (it, uh, φh) in optimiser
     write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)

scripts/MPI/3d_thermal_compliance_ALM.jl

@@ -3,12 +3,12 @@ using Gridap, Gridap.MultiField, GridapDistributed, GridapPETSc, GridapSolvers,
 
 using GridapSolvers: NewtonSolver
 
-global elx = parse(Int,ARGS[1])
-global ely = parse(Int,ARGS[2])
-global elz = parse(Int,ARGS[3])
-global Px = parse(Int,ARGS[4])
-global Py = parse(Int,ARGS[5])
-global Pz = parse(Int,ARGS[6])
+# global elx = parse(Int,ARGS[1])
+# global ely = parse(Int,ARGS[2])
+# global elz = parse(Int,ARGS[3])
+# global Px = parse(Int,ARGS[4])
+# global Py = parse(Int,ARGS[5])
+# global Pz = parse(Int,ARGS[6])
 
 """
   (MPI) Minimum thermal compliance with augmented Lagrangian method in 3D.
@@ -73,7 +73,7 @@ function main(mesh_partition,distribute,el_size,coef,step)
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
   a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  l(v,φ,dΩ,dΓ_N) = ∫(10v)dΓ_N
 
   ## Optimisation functionals
   J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
@@ -119,19 +119,19 @@ function main(mesh_partition,distribute,el_size,coef,step)
 end
 
 with_mpi() do distribute
-  mesh_partition = (Px,Py,Pz)
-  el_size = (elx,ely,elz)
+  mesh_partition = (5,5,5)#(Px,Py,Pz)
+  el_size = (150,150,150)#(elx,ely,elz)
   all_solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true 
     -ksp_converged_reason -ksp_rtol 1.0e-12"
 
   GridapPETSc.with(args=split(all_solver_options)) do
     main(mesh_partition,distribute,el_size,5,10)
-    main(mesh_partition,distribute,el_size,2,10)
+    # main(mesh_partition,distribute,el_size,2,10)
 
-    main(mesh_partition,distribute,el_size,5,5)
-    main(mesh_partition,distribute,el_size,2,5)
+    # main(mesh_partition,distribute,el_size,5,5)
+    # main(mesh_partition,distribute,el_size,2,5)
 
-    main(mesh_partition,distribute,el_size,5,3)
-    main(mesh_partition,distribute,el_size,2,3)
+    # main(mesh_partition,distribute,el_size,5,3)
+    # main(mesh_partition,distribute,el_size,2,3)
   end
 end;

scripts/MPI/inverse_homenisation_ALM.jl

@@ -44,7 +44,7 @@ function main(mesh_partition,distribute,el_size)
   reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
   reffe_scalar = ReferenceFE(lagrangian,Float64,order)
   V = TestFESpace(model,reffe;dirichlet_tags=["origin"])
-  U = TrialFESpace(_V,VectorValue(0.0,0.0))
+  U = TrialFESpace(V,VectorValue(0.0,0.0))
   V_reg = V_φ = TestFESpace(model,reffe_scalar)
   U_reg = TrialFESpace(V_reg)
 
@@ -64,12 +64,10 @@ function main(mesh_partition,distribute,el_size)
   l = [(v,φ,dΩ) -> ∫(-(I ∘ φ) * C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
 
   ## Optimisation functionals
-  _C(C,ε_p,ε_q) = C ⊙ ε_p ⊙ ε_q;
-  _K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],εᴹ[2]))/4
-  _v_K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],ε(u1)+εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],ε(u2)+εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],ε(u2)+εᴹ[2]))/4    
-
-  J(u,φ,dΩ) = ∫(-(I ∘ φ)*_K(C,u,εᴹ))dΩ
-  dJ(q,u,φ,dΩ) = ∫(_v_K(C,u,εᴹ)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ))
+  dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ))
   Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
   dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
@@ -88,7 +86,7 @@ function main(mesh_partition,distribute,el_size)
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
     ls = solver,adjoint_ls = solver
   )
-  pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dJ=dJ,analytic_dC=[dVol])
+  pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
   α = α_coeff*maximum(el_Δ)

scripts/MPI/jobs/3d_thermal_compliance_ALM.pbs

@@ -4,19 +4,23 @@
 #PBS -N 3d_thermal_compliance_ALM
 
 #PBS -l select=125:ncpus=1:mpiprocs=1:ompthreads=1:mem=8GB:cputype=6140
-#PBS -l walltime=200:00:00
+#PBS -l walltime=23:00:00
 #PBS -j oe
 
 source $HOME/hpc-environments-main/lyra/load-ompi.sh
 PROJECT_DIR=$HOME/LevelSetTopOpt/
 
 julia --project=$PROJECT_DIR -e "using Pkg; Pkg.precompile()"
 
-mpirun --hostfile $PBS_NODEFILE julia --project=$PROJECT_DIR --check-bounds no -O3 --compiled-modules=no \
-    $PROJECT_DIR/scripts/MPI/3d_thermal_compliance_ALM.jl \
-    150 \
-    150 \
-    150 \
-    5 \
-    5 \
-    5
+mpirun --hostfile $PBS_NODEFILE \
+ julia --project=$PROJECT_DIR --check-bounds no -O3 --compiled-modules=no \
+  $PROJECT_DIR/scripts/MPI/3d_thermal_compliance_ALM.jl
+
+# mpirun --hostfile $PBS_NODEFILE julia --project=$PROJECT_DIR --check-bounds no -O3 --compiled-modules=no \
+#     $PROJECT_DIR/scripts/MPI/3d_thermal_compliance_ALM.jl \
+#     150 \
+#     150 \
+#     150 \
+#     5 \
+#     5 \
+#     5

scripts/Serial/inverse_homenisation_ALM.jl

@@ -62,12 +62,10 @@ function main()
   l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
 
   ## Optimisation functionals
-  _C(C,ε_p,ε_q) = C ⊙ ε_p ⊙ ε_q
-  _K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],εᴹ[2]))/4
-  _v_K(C,(u1,u2,u3),εᴹ) = (_C(C,ε(u1)+εᴹ[1],ε(u1)+εᴹ[1]) + _C(C,ε(u2)+εᴹ[2],ε(u2)+εᴹ[2]) + 2*_C(C,ε(u1)+εᴹ[1],ε(u2)+εᴹ[2]))/4    
-
-  J(u,φ,dΩ) = ∫(-(I ∘ φ)*_K(C,u,εᴹ))dΩ
-  dJ(q,u,φ,dΩ) = ∫(_v_K(C,u,εᴹ)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ))
+  dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ))
   Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
   dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
@@ -76,7 +74,7 @@ function main()
 
   ## Setup solver and FE operators
   state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ)
-  pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dJ=dJ,analytic_dC=[dVol])
+  pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
   α = α_coeff*maximum(el_Δ)

scripts/_dev/bug_issue_46/thermal_MPI.jl

@@ -0,0 +1,91 @@
+using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays, 
+  SparseMatricesCSR
+
+function main(mesh_partition,distribute,use_l::Bool)
+  ranks = distribute(LinearIndices((prod(mesh_partition),)))
+  # FE parameters
+  order = 1                                       # Finite element order
+  xmax = ymax = 1.0                               # Domain size
+  dom = (0,xmax,0,ymax)                           # Bounding domain
+  el_size = (200,200)                             # Mesh partition size
+  prop_Γ_N = 0.2                                  # Γ_N size parameter
+  prop_Γ_D = 0.2                                  # Γ_D size parameter
+  f_Γ_N(x) = (x[1] ≈ xmax &&                      # Γ_N indicator function
+    ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + eps())
+  f_Γ_D(x) = (x[1] ≈ 0.0 &&                       # Γ_D indicator function
+    (x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps())) # @$\lvert\lvert$@
+  # FD parameters
+  γ = 0.1                                         # HJ eqn time step coeff
+  γ_reinit = 0.5                                  # Reinit. eqn time step coeff
+  max_steps = floor(Int,minimum(el_size)/10)      # Max steps for advection
+  tol = 1/(5order^2)/minimum(el_size)             # Advection tolerance
+  # Problem parameters
+  κ = 1                                           # Diffusivity
+  g = 1                                           # Heat flow in
+  vf = 0.4                                        # Volume fraction constraint
+  lsf_func = initial_lsf(4,0.2)                   # Initial level set function
+  iter_mod = 10                                   # VTK Output modulo
+  path = "./results/AD_tut1_MPI_lobj_$use_l/"              # Output path
+  i_am_main(ranks) && mkpath(path)                # Create path
+  # Model
+  model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size);
+  update_labels!(1,model,f_Γ_D,"Gamma_D")
+  update_labels!(2,model,f_Γ_N,"Gamma_N")
+  # Triangulation and measures
+  Ω = Triangulation(model)
+  Γ_N = BoundaryTriangulation(model,tags="Gamma_N")
+  dΩ = Measure(Ω,2*order)
+  dΓ_N = Measure(Γ_N,2*order)
+  # Spaces
+  reffe = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"])
+  U = TrialFESpace(V,0.0)
+  V_φ = TestFESpace(model,reffe)
+  V_reg = TestFESpace(model,reffe;dirichlet_tags=["Gamma_N"])
+  U_reg = TrialFESpace(V_reg,0)
+  # Level set and interpolator
+  φh = interpolate(lsf_func,V_φ)
+  interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model)))
+  I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
+  # Weak formulation
+  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ,dΩ,dΓ_N) = ∫(g*v)dΓ_N
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  # Objective and constraints
+  J(u,φ,dΩ,dΓ_N) = if use_l
+    ∫(g*u)dΓ_N #+ ∫(0)dΩ
+  else
+    ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  end
+  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  vol_D = sum(∫(1)dΩ)
+  C(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dC(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  pcfs = PDEConstrainedFunctionals(J,[C],state_map)
+  # Velocity extension
+  α = 4*maximum(get_el_Δ(model))
+  a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ
+  vel_ext = VelocityExtension(a_hilb, U_reg, V_reg)
+  # Finite difference scheme
+  scheme = FirstOrderStencil(length(el_size),Float64)
+  stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
+  # Optimiser
+  optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
+    γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
+  # Solve
+  for (it,uh,φh) in optimiser
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && (writevtk(Ω,path*"out$it",cellfields=data);GC.gc())
+    write_history(path*"/history.txt",get_history(optimiser);ranks)
+  end
+  # Final structure
+  it = get_history(optimiser).niter; uh = get_state(pcfs)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,
+    "H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+end
+
+with_mpi() do distribute
+  mesh_partition = (2,2)
+  main(mesh_partition,distribute,false)
+  main(mesh_partition,distribute,true) # <- this should crash
+end