Towards distributed

scripts/Embedded/Bugs/FIXED_change_domain_bug.jl

@@ -25,7 +25,7 @@ dΓ = Measure(Γ,2order)
 
 
 _Γ = EmbeddedBoundary(cutgeo)
-Γ = DifferentiableTriangulation(_Γ)
+Γ = DifferentiableTriangulation(_Γ,V_φ)
 dΓ = Measure(Γ,2order)
 ∫(vh)dΓ
 ∫(2vh)dΓ

scripts/Embedded/Bugs/basis_isbitstype_bug.jl

@@ -41,7 +41,7 @@ cutgeo = cut(model,geo)
 # A.1) Volume integral
 
 Ω = Triangulation(cutgeo,PHYSICAL_IN)
-Ω_AD = DifferentiableTriangulation(Ω)
+Ω_AD = DifferentiableTriangulation(Ω,V_φ)
 dΩ = Measure(Ω_AD,2*order)
 
 Γ = EmbeddedBoundary(cutgeo)

scripts/Embedded/Bugs/interface_change_domain_bug.jl

@@ -0,0 +1,45 @@
+using Test
+
+using GridapTopOpt
+using Gridap
+
+using GridapDistributed, PartitionedArrays
+
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Adaptivity
+
+using GridapTopOpt: get_subfacet_normal_vector, get_ghost_normal_vector
+using GridapTopOpt: get_conormal_vector, get_tangent_vector
+
+using GridapDistributed: DistributedTriangulation, DistributedDomainContribution
+
+order = 1
+n = 16
+
+parts = (2,2)
+ranks = DebugArray(LinearIndices((prod(parts),)))
+
+# _model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(n,n))
+_model = CartesianDiscreteModel((0,1,0,1),(n,n))
+
+base_model = UnstructuredDiscreteModel(_model)
+ref_model = refine(base_model, refinement_method = "barycentric")
+model = Gridap.Adaptivity.get_model(ref_model)
+
+reffe = ReferenceFE(lagrangian,Float64,order)
+V_φ = TestFESpace(model,reffe)
+
+φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.5223,V_φ) # Circle
+fh = interpolate(x->cos(x[1]*x[2]),V_φ)
+
+geo = DiscreteGeometry(φh,model)
+cutgeo = cut(model,geo)
+
+Γ = EmbeddedBoundary(cutgeo)
+Γ_AD = DifferentiableTriangulation(Γ,V_φ)
+dΓ_AD = Measure(Γ_AD,2*order)
+
+J_int2(φ) = ∫(g(fh))dΓ_AD
+dJ_int_AD2 = gradient(J_int2,φh)
+dJ_int_AD_vec2 = assemble_vector(dJ_int_AD2,V_φ)

scripts/Embedded/Examples/2d_compliance_L.jl

@@ -56,11 +56,11 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
   n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -70,7 +70,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 function lame_parameters(E,ν)
@@ -87,7 +87,7 @@ E = 1.0
 γg = 0.1
 
 g = VectorValue(0,-1)
-a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin + 
+a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin +
   # ∫((γg*h)*jump(nΛ_D⋅∇(v)) ⋅ jump(nΛ_D⋅∇(u)))dΛ_D + # <- this currently breaks due to `_compute_cell_ids` function
   ∫((γg*h^3)*jump(Ωs.n_Γg⋅∇(v)) ⋅ jump(Ωs.n_Γg⋅∇(u)))Ωs.dΓg
 l(v,φ) = ∫(v⋅g)dΓ_N
@@ -106,12 +106,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,VectorValue(0.0,0.0))
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method

scripts/Embedded/Examples/2d_thermal.jl

@@ -44,11 +44,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -58,7 +57,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -78,12 +77,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method

scripts/Embedded/Examples/2d_thermal_AgFEM.jl

@@ -44,11 +44,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -58,7 +57,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -81,12 +80,12 @@ state_collection = EmbeddedCollection(model,φh) do cutgeo,_
   V = AgFEMSpace(Vstd,aggregates)
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method

scripts/Embedded/Examples/2d_thermal_L.jl

@@ -52,11 +52,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -66,7 +65,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -86,12 +85,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method

scripts/Embedded/Examples/2d_thermal_gmsh_example.jl

@@ -52,11 +52,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.3,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -66,7 +65,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -86,12 +85,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method

scripts/Embedded/Examples/2d_thermal_weird_background.jl

@@ -52,11 +52,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.3,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -66,7 +65,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -86,12 +85,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method

scripts/Embedded/Examples/FCM_2d_compliance_L.jl

@@ -55,9 +55,8 @@ V_φ = TestFESpace(model,reffe_scalar)
   Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
   Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Ωout  => Ωout,
@@ -69,7 +68,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 function lame_parameters(E,ν)
@@ -84,7 +83,7 @@ E = 1.0
 g = VectorValue(0,-0.1)
 
 const ϵ = (λ + μ)*1e-3
-a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin + 
+a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin +
   # ∫(ϵ*(ε(v) ⊙ (σ ∘ ε(u))))Ωs.dΩout # Ersatz
   ∫(ϵ * ∇(v) ⊙ ∇(u))Ωs.dΩout # Pure diffusion
 l(v,φ) = ∫(v⋅g)dΓ_N

scripts/Embedded/Examples/FCM_2d_thermal.jl

@@ -47,9 +47,8 @@ V_φ = TestFESpace(model,reffe_scalar)
   Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
   Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Ωout  => Ωout,

scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl

@@ -4,7 +4,10 @@ using GridapEmbedded, GridapEmbedded.LevelSetCutters
 
 using GridapTopOpt: StateParamIntegrandWithMeasure
 
-using GridapDistributed, PartitionedArrays, GridapPETSc
+using GridapDistributed, GridapPETSc, GridapSolvers,
+  PartitionedArrays, GridapTopOpt, SparseMatricesCSR
+
+using GridapSolvers: NewtonSolver
 
 function main(mesh_partition,distribute,el_size,path)
   ranks = distribute(LinearIndices((prod(mesh_partition),)))
@@ -53,9 +56,8 @@ function main(mesh_partition,distribute,el_size,path)
     Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
     Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
     Γg = GhostSkeleton(cutgeo)
-    n_Γg = get_normal_vector(Γg)
     Ωact = Triangulation(cutgeo,ACTIVE)
-    (; 
+    (;
       :Ωin  => Ωin,
       :dΩin => Measure(Ωin,2*order),
       :Ωout  => Ωout,
@@ -95,8 +97,13 @@ function main(mesh_partition,distribute,el_size,path)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dC=(dVol,))
 
   ## Evolution Method
-  evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps)
-  reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0))
+  evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps,
+    ode_ls = solver,
+    assembler=SparseMatrixAssembler(Tm,Tv,V_φ,V_φ))
+  reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;
+    stabilisation_method=ArtificialViscosity(3.0),
+    assembler=SparseMatrixAssembler(Tm,Tv,V_φ,V_φ),
+    nls = NewtonSolver(LUSolver();maxiter=50,rtol=1.e-14,verbose=i_am_main(ranks)))
   ls_evo = UnfittedFEEvolution(evo,reinit)
   reinit!(ls_evo,φh)
 
@@ -118,9 +125,6 @@ function main(mesh_partition,distribute,el_size,path)
   optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;debug=true,
     γ,verbose=i_am_main(ranks),constraint_names=[:Vol],converged)
   for (it,uh,φh,state) in optimiser
-    x_φ = get_free_dof_values(φh)
-    idx = findall(isapprox(0.0;atol=10^-10),x_φ)
-    !isempty(idx) && i_am_main(ranks) && @warn "Boundary intersects nodes!"
     if iszero(it % iter_mod)
       writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"velh"=>FEFunction(V_φ,state.vel)])
       writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])