Merge pull request #274 from TzuYaoHuang/fix/DotInPeriodicBC

Fix Gram-Schmidt parameter in PCG under periodic BC

Author: weymouth

src/Poisson.jl

@@ -129,7 +129,7 @@ function pcg!(p::Poisson{T};it=6) where T
     for i in 1:it
         perBC!(ϵ,p.perdir)
         @inside z[I] = mult(I,p.L,p.D,ϵ)
-        alpha = rho/(z⋅ϵ)
+        alpha = rho/perdot(z,ϵ,p.perdir)
         (abs(alpha)<1e-2 || abs(alpha)>1e2) && return # alpha should be O(1)
         @loop (x[I] += alpha*ϵ[I];
                r[I] -= alpha*z[I]) over I ∈ inside(x)

src/util.jl

@@ -257,6 +257,14 @@ perBC!(a, perdir, N = size(a)) = for j ∈ perdir
     @loop a[I] = a[CIj(j,I,N[j]-1)] over I ∈ slice(N,1,j)
     @loop a[I] = a[CIj(j,I,2)] over I ∈ slice(N,N[j],j)
 end
+using LinearAlgebra: ⋅
+"""
+    perdot(a,b,perdir)
+
+Apply dot product to the inner cells of two _scalar_ fields, assuming zero values in ghost cell when using Neumann BC.
+"""
+perdot(a,b,::Tuple{}) = a⋅b
+perdot(a,b,perdir,R=inside(a)) = @view(a[R])⋅@view(b[R])
 """
     interp(x::SVector, arr::AbstractArray)
 

test/maintests.jl

@@ -94,6 +94,14 @@ backend != "KernelAbstractions" && throw(ArgumentError("SIMD backend not allowed
        @test GPUArrays.@allowscalar WaterLily.interp(SVector(3.5,3),b) ≈ 3.5
        @test GPUArrays.@allowscalar eltype(WaterLily.interp(SVector(3.5,3),b))==Float64
        @test_throws MethodError GPUArrays.@allowscalar WaterLily.interp(SVector(2.5f0,1.f0),b)
+
+        # test on perdot
+        σ1 = rand(Ng...) |> f # scalar
+        σ2 = rand(Ng...) |> f # another scalar 
+        # use ≈ instead of == as summation in different order might result in slight difference in floating point expressions
+        @test GPUArrays.@allowscalar WaterLily.perdot(σ1,σ2,())    ≈ sum(σ1[I]*σ2[I] for I∈CartesianIndices(σ1))
+        @test GPUArrays.@allowscalar WaterLily.perdot(σ1,σ2,(1,))  ≈ sum(σ1[I]*σ2[I] for I∈inside(σ1))
+        @test GPUArrays.@allowscalar WaterLily.perdot(σ1,σ2,(1,2)) ≈ sum(σ1[I]*σ2[I] for I∈inside(σ1))
     end
 end