From d06626b89543ad7b256166793894941671223042 Mon Sep 17 00:00:00 2001
From: Pere Antoni Martorell <pmartorell@cimne.upc.edu>
Date: Tue, 30 Apr 2024 09:25:06 +0200
Subject: [PATCH] [distributed] add distributed poisson tests

---
 test/DistributedTests/CutterTests.jl          | 10 +-
 test/DistributedTests/Poisson.jl              | 98 +++++++++++++++++++
 test/DistributedTests/mpi/runtests_body.jl    |  6 ++
 .../sequential/CutterTests.jl                 |  2 +
 .../sequential/PoissonTests.jl                | 16 +++
 5 files changed, 127 insertions(+), 5 deletions(-)
 create mode 100644 test/DistributedTests/Poisson.jl
 create mode 100644 test/DistributedTests/sequential/PoissonTests.jl

diff --git a/test/DistributedTests/CutterTests.jl b/test/DistributedTests/CutterTests.jl
index e1c86e03..8b22ec7d 100644
--- a/test/DistributedTests/CutterTests.jl
+++ b/test/DistributedTests/CutterTests.jl
@@ -68,17 +68,17 @@ function main(distribute;
     writevtk(Γ,"Γ")
   end
 
-  @test vol_in + vol_out - vol_box < 1e-10
+  @test vol_in + vol_out - vol_box < 1e-9
   ref_vol,ref_surf = reference_volume_and_surface(geoname)
   if ref_vol > 0
-    @test abs(vol_in - ref_vol) < 1e-10
-    @test abs(surf - ref_surf) < 1e-10
+    @test abs(vol_in - ref_vol) < 1e-9
+    @test abs(surf - ref_surf) < 1e-9
   end
 end
 
 function reference_volume_and_surface(geoname)
-  volumes = Dict("cube"=>1.0,)
-  surfaces = Dict("cube"=>6.0,)
+  volumes = Dict("cube"=>1.0,"Bunny-LowPoly"=>273280.03374196636)
+  surfaces = Dict("cube"=>6.0,"Bunny-LowPoly"=>29490.7154966073)
   if haskey(volumes,geoname)
     volumes[geoname],surfaces[geoname]
   else
diff --git a/test/DistributedTests/Poisson.jl b/test/DistributedTests/Poisson.jl
new file mode 100644
index 00000000..82805f7b
--- /dev/null
+++ b/test/DistributedTests/Poisson.jl
@@ -0,0 +1,98 @@
+module DistributedPoissonTests
+
+using Gridap
+using STLCutters
+using GridapDistributed
+using PartitionedArrays
+using GridapEmbedded
+using Test
+
+function main(distribute;
+  np = (1,1,1),
+  nc = (2,2,2),
+  geoname = "cube",
+  δ = 0.2,
+  tolfactor=10^4,
+  simplex=false,
+  vtk = false,
+  verbose = false)
+
+  # Manufactured solution
+  u(x) = x[1] + x[2] - x[3]
+  f(x) = - Δ(u)(x)
+  ud(x) = u(x)
+
+  filename = joinpath(@__DIR__,"..","data","$geoname.stl")
+  if !isfile(filename)
+    filename = download_thingi10k(id;path="")
+  end
+
+  ranks = distribute(LinearIndices((prod(np),)))
+  geo = STLGeometry(filename)
+
+  pmin,pmax = get_bounding_box(geo)
+  diagonal = pmax-pmin
+  pmin = pmin - diagonal*δ
+  pmax = pmax + diagonal*δ
+  bgmodel = CartesianDiscreteModel(ranks,np,pmin,pmax,nc)
+  if simplex
+    bgmodel = simplexify(bgmodel,positive=true)
+  end
+
+  cutter = STLCutter(;tolfactor)
+  cutgeo = cut(cutter,bgmodel,geo)
+
+  Ω = Triangulation(cutgeo,PHYSICAL_IN)
+  Γ = EmbeddedBoundary(cutgeo)
+  n_Γ = get_normal_vector(Γ)
+  degree = 2
+  dΩ = Measure(Ω,degree)
+  dΓ = Measure(Γ,degree)
+
+  # Setup FESpace
+  order = 1
+  Ω_act = Triangulation(cutgeo,ACTIVE)
+  V = TestFESpace(Ω_act,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
+  U = TrialFESpace(V)
+
+  # Weak form
+  γ = 10.0
+  h = (pmax - pmin)[1] / nc[1]
+
+  a(u,v) =
+    ∫( ∇(v)⋅∇(u) ) * dΩ +
+    ∫( (γ/h)*v*u  - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u ) * dΓ
+
+  l(v) =
+    ∫( v*f ) * dΩ +
+    ∫( (γ/h)*v*ud - (n_Γ⋅∇(v))*ud ) * dΓ
+
+  # FE problem
+  @time op = AffineFEOperator(a,l,U,V)
+  @time uh = solve(op)
+
+  e = u - uh
+
+  # Postprocess
+  l2(u) = sqrt(sum( ∫( u*u )*dΩ ))
+  h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ ))
+
+  el2 = l2(e)
+  eh1 = h1(e)
+  ul2 = l2(uh)
+  uh1 = h1(uh)
+
+  if vtk
+    writevtk(Ω,"results",cellfields=["uh"=>uh])
+  end
+
+  if verbose
+    println("L2 error: ", el2/ul2)
+    println("H1 error: ", eh1/uh1)
+  end
+
+  @test el2/ul2 < 1.e-8
+  @test eh1/uh1 < 1.e-7
+end
+
+end # module
diff --git a/test/DistributedTests/mpi/runtests_body.jl b/test/DistributedTests/mpi/runtests_body.jl
index 754d2fcc..f2c4af18 100644
--- a/test/DistributedTests/mpi/runtests_body.jl
+++ b/test/DistributedTests/mpi/runtests_body.jl
@@ -5,6 +5,7 @@ const PArrays = PartitionedArrays
 using MPI
 
 include("../CutterTests.jl")
+include("../Poisson.jl")
 
 if ! MPI.Initialized()
   MPI.Init()
@@ -22,6 +23,11 @@ function all_tests(distribute,parts)
   DistributedCutterTests.main(distribute,np=parts,nc=(8,8,8),simplex=true)
   PArrays.toc!(t,"MPIDistributedCutter")
 
+  DistributedPoissonTests.main(distribute,np=parts,nc=(4,4,4))
+  DistributedPoissonTests.main(distribute,np=parts,nc=(8,8,8))
+  DistributedPoissonTests.main(distribute,np=(2,2,2),nc=(8,8,8),geoname="Bunny-LowPoly")
+  PArrays.toc!(t,"MPIDistributedPoisson")
+
   display(t)
 end
 
diff --git a/test/DistributedTests/sequential/CutterTests.jl b/test/DistributedTests/sequential/CutterTests.jl
index 55b15655..4f48e069 100644
--- a/test/DistributedTests/sequential/CutterTests.jl
+++ b/test/DistributedTests/sequential/CutterTests.jl
@@ -14,10 +14,12 @@ with_debug() do distribute
   main(distribute,np=(2,2,2),nc=(8,8,8))
   main(distribute,np=(2,2,2),nc=(4,4,4),simplex=true)
   main(distribute,np=(2,2,2),nc=(8,8,8),simplex=true)
+  main(distribute,np=(2,2,2),nc=(8,8,8),geoname="Bunny-LowPoly")
 
   # Distributed with propagation
   main(distribute,np=(3,3,3),nc=(9,9,9))
   main(distribute,np=(3,3,3),nc=(9,9,9),simplex=true)
+  main(distribute,np=(3,3,3),nc=(9,9,9),geoname="Bunny-LowPoly")
 
   # Untouched subdomain interiors
   main(distribute,np=(3,3,3),nc=(9,9,9),δ=0.5)
diff --git a/test/DistributedTests/sequential/PoissonTests.jl b/test/DistributedTests/sequential/PoissonTests.jl
new file mode 100644
index 00000000..ee41c8dc
--- /dev/null
+++ b/test/DistributedTests/sequential/PoissonTests.jl
@@ -0,0 +1,16 @@
+module PoissonTestsSeq
+using PartitionedArrays
+include("../Poisson.jl")
+main = DistributedPoissonTests.main
+
+with_debug() do distribute
+
+  # Degenerated Case
+  main(distribute,np=(1,1,1),nc=(4,4,4))
+
+  # Distributed without propagation
+  main(distribute,np=(2,2,2),nc=(4,4,4))
+  main(distribute,np=(2,2,2),nc=(8,8,8))
+
+end
+end