From 0682fc913badb552ad924f38d9fedd0873dc3c94 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> Date: Thu, 4 Jul 2024 22:42:31 +0200 Subject: [PATCH] Add missing 3D no-slip boundary condition (#187) * add missing 3D no_slip * 3D BCs tests * 3D no_slip * docs * format * fix CI * fix no_slip in 3D * remove commented code * remove commented code * remove file * fix 3d velocity interpolation * small bug * format * remove old lines of code * format * fix CI? * format * fix * fix again * format a bit * fix no slip2D --- docs/make.jl | 1 + docs/src/man/boundary_conditions.md | 36 +++++++++ src/Interpolations.jl | 24 ++++-- src/boundaryconditions/BoundaryConditions.jl | 7 +- src/boundaryconditions/no_slip.jl | 77 ++++++++++++------ src/stokes/Stokes3D.jl | 17 ++-- test/test_boundary_conditions3D.jl | 84 ++++++++++++++++++++ 7 files changed, 208 insertions(+), 38 deletions(-) create mode 100644 docs/src/man/boundary_conditions.md create mode 100644 test/test_boundary_conditions3D.jl diff --git a/docs/make.jl b/docs/make.jl index e34f462f..0824923f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,6 +16,7 @@ makedocs(; "Installation" => "man/installation.md", "Backend" => "man/backend.md", "Equations" => "man/equations.md", + "Boundary conditions" => "man/boundary_conditions.md", "Advection" => "man/advection.md", ], "Examples" => Any[ diff --git a/docs/src/man/boundary_conditions.md b/docs/src/man/boundary_conditions.md new file mode 100644 index 00000000..98665708 --- /dev/null +++ b/docs/src/man/boundary_conditions.md @@ -0,0 +1,36 @@ +# Flow boundary conditions + +Supported boundary conditions: + +1. Free slip + + $\frac{\partial u_i}{\partial x_i} = 0$ at the boundary $\Gamma$ + +2. No slip + + $u_i = 0$ at the boundary $\Gamma$ + +3. Free surface + + $\sigma_z = 0 \rightarrow \tau_z = P$ at the top boundary + +## Defining the boundary contions +Information regarding flow boundary conditions is defined in the `FlowBoundaryConditions` object. They can be switched on and off by setting them as `true` or `false` at the appropriate boundaries. Valid boundary names are `left` and `right`, `top` and `bot`, and for the 3D case, `front` and `back`. + +For example, if we want to have free free-slip in every single boundary in a 2D simulation, we need to instantiate `FlowBoundaryConditions` as: +```julia +bcs = FlowBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false), + free_slip = (left=true, right=true, top=true, bot=true), + free_surface = false +) +``` + +The equivalent for the 3D case would be: +```julia +bcs = FlowBoundaryConditions(; + no_slip = (left=false, right=false, top=false, bot=false, front=false, back=false), + free_slip = (left=true, right=true, top=true, bot=true, front=true, back=true), + free_surface = false +) +``` \ No newline at end of file diff --git a/src/Interpolations.jl b/src/Interpolations.jl index f6b5a102..c681e71b 100644 --- a/src/Interpolations.jl +++ b/src/Interpolations.jl @@ -258,21 +258,31 @@ onto the pre-allocated `Vx_d`, `Vy_d`, `Vz_d` 3D arrays located at the grid vert function velocity2vertex!(Vx_v, Vy_v, Vz_v, Vx, Vy, Vz) # infer size of grid nx, ny, nz = size(Vx) + n = max(nx, ny, nz) nv_x, nv_y, nv_z = nx - 1, ny - 2, nz - 2 # interpolate to cell vertices - @parallel (@idx nv_x nv_y nv_z) _velocity2vertex!(Vx_v, Vy_v, Vz_v, Vx, Vy, Vz) + @parallel (@idx n, n, n) _velocity2vertex!(Vx_v, Vy_v, Vz_v, Vx, Vy, Vz) return nothing end @parallel_indices (i, j, k) function _velocity2vertex!(Vx_v, Vy_v, Vz_v, Vx, Vy, Vz) @inbounds begin - Vx_v[i, j, k] = - 0.25 * (Vx[i, j, k] + Vx[i, j + 1, k] + Vx[i, j, k + 1] + Vx[i, j + 1, k + 1]) - Vy_v[i, j, k] = - 0.25 * (Vy[i, j, k] + Vy[i + 1, j, k] + Vy[i, j, k + 1] + Vy[i + 1, j, k + 1]) - Vz_v[i, j, k] = - 0.25 * (Vz[i, j, k] + Vz[i, j + 1, k] + Vz[i + 1, j, k] + Vz[i + 1, j + 1, k]) + if all((i, j, k) .≤ size(Vx)) + Vx_v[i, j, k] = + 0.25 * + (Vx[i, j, k] + Vx[i, j + 1, k] + Vx[i, j, k + 1] + Vx[i, j + 1, k + 1]) + end + if all((i, j, k) .≤ size(Vy)) + Vy_v[i, j, k] = + 0.25 * + (Vy[i, j, k] + Vy[i + 1, j, k] + Vy[i, j, k + 1] + Vy[i + 1, j, k + 1]) + end + if all((i, j, k) .≤ size(Vz)) + Vz_v[i, j, k] = + 0.25 * + (Vz[i, j, k] + Vz[i, j + 1, k] + Vz[i + 1, j, k] + Vz[i + 1, j + 1, k]) + end end return nothing end diff --git a/src/boundaryconditions/BoundaryConditions.jl b/src/boundaryconditions/BoundaryConditions.jl index 81968a6a..011039d3 100644 --- a/src/boundaryconditions/BoundaryConditions.jl +++ b/src/boundaryconditions/BoundaryConditions.jl @@ -53,7 +53,12 @@ end function _flow_bcs!(bcs, V) n = bc_index(V) # no slip boundary conditions - do_bc(bcs.no_slip) && (@parallel (@idx n) no_slip!(V..., bcs.no_slip)) + # do_bc(bcs.no_slip) && (@parallel (@idx n) no_slip!(V..., bcs.no_slip)) + if do_bc(bcs.no_slip) + # @parallel (@idx n) no_slip1!(V..., bcs.no_slip) + # @parallel (@idx n) no_slip2!(V..., bcs.no_slip) + no_slip!(V..., bcs.no_slip) + end # free slip boundary conditions do_bc(bcs.free_slip) && (@parallel (@idx n) free_slip!(V..., bcs.free_slip)) diff --git a/src/boundaryconditions/no_slip.jl b/src/boundaryconditions/no_slip.jl index 2a677277..7295d797 100644 --- a/src/boundaryconditions/no_slip.jl +++ b/src/boundaryconditions/no_slip.jl @@ -1,25 +1,54 @@ -@parallel_indices (i) function no_slip!(Ax, Ay, bc) - @inbounds begin - if bc.left - (i ≤ size(Ax, 2)) && (Ax[1, i] = 0.0) - (1 < i < size(Ay, 2)) && (Ay[1, i] = -Ay[2, i]) - end - if bc.right - (i ≤ size(Ax, 2)) && (Ax[end, i] = 0.0) - (1 < i < size(Ay, 2)) && (Ay[end, i] = -Ay[end - 1, i]) - end - if bc.bot - (i ≤ size(Ay, 1)) && (Ay[i, 1] = 0.0) - (1 < i < size(Ax, 1)) && (Ax[i, 1] = -Ax[i, 2]) - end - if bc.top - (i ≤ size(Ay, 1)) && (Ay[i, end] = 0.0) - (1 < i < size(Ax, 1)) && (Ax[i, end] = -Ax[i, end - 1]) - end - # corners - # bc.bot && (Ax[1, 1] = 0.0; Ax[1, 1] = 0.0) - # bc.left && bc.bot && (Ax[1, 1] = 0.0) - # bc.right && bc.top && (Ay[end, end] = 0.0) - end - return nothing +@views function no_slip!(Ax, Ay, bc) + if bc.left + Ax[1, :] .= 0 + Ay[1, :] .= -Ay[2, :] + end + if bc.right + Ax[end, :] .= 0 + Ay[end, :] .= -Ay[end - 1, :] + end + if bc.bot + Ax[:, 1] .= -Ax[:, 2] + Ay[:, 1] .= 0 + end + if bc.top + Ax[:, end] .= -Ax[:, end - 1] + Ay[:, end] .= 0 + end +end + +@views function no_slip!(Ax, Ay, Az, bc) + if bc.left + Ax[1, :, :] .= 0 + Ay[1, :, :] .= -Ay[2, :, :] + Az[1, :, :] .= -Az[2, :, :] + end + if bc.right + Ax[end, :, :] .= 0 + Ay[end, :, :] .= -Ay[end - 1, :, :] + Az[end, :, :] .= -Az[end - 1, :, :] + end + + if bc.front + Ax[:, 1, :] .= -Ax[:, 2, :] + Ay[:, 1, :] .= 0 + Az[:, 1, :] .= -Az[:, 2, :] + end + + if bc.back + Ax[:, end, :] .= -Ax[:, end - 1, :] + Ay[:, end, :] .= 0 + Az[:, end, :] .= -Az[:, end - 1, :] + end + + if bc.bot + Ax[:, :, 1] .= -Ax[:, :, 2] + Ay[:, :, 1] .= -Ay[:, :, 2] + Az[:, :, 1] .= 0 + end + if bc.top + Ax[:, :, end] .= -Ax[:, :, end - 1] + Ay[:, :, end] .= -Ay[:, :, end - 1] + Az[:, :, end] .= 0 + end end diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index b4ff2ee3..c8b1975e 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -292,10 +292,15 @@ function _solve!( iter += 1 if iter % nout == 0 && iter > 1 cont += 1 - push!(norm_Rx, maximum_mpi(abs.(stokes.R.Rx))) - push!(norm_Ry, maximum_mpi(abs.(stokes.R.Ry))) - push!(norm_Rz, maximum_mpi(abs.(stokes.R.Rz))) - push!(norm_∇V, maximum_mpi(abs.(stokes.R.RP))) + # push!(norm_Rx, maximum_mpi(abs.(stokes.R.Rx))) + # push!(norm_Ry, maximum_mpi(abs.(stokes.R.Ry))) + # push!(norm_Rz, maximum_mpi(abs.(stokes.R.Rz))) + # push!(norm_∇V, maximum_mpi(abs.(stokes.R.RP))) + push!(norm_Rx, norm_mpi(stokes.R.Rx) / length(stokes.R.Rx)) + push!(norm_Ry, norm_mpi(stokes.R.Ry) / length(stokes.R.Ry)) + push!(norm_Rz, norm_mpi(stokes.R.Rz) / length(stokes.R.Rz)) + push!(norm_∇V, norm_mpi(stokes.R.RP) / length(stokes.R.RP)) + err = max(norm_Rx[cont], norm_Ry[cont], norm_Rz[cont], norm_∇V[cont]) push!(err_evo1, err) push!(err_evo2, iter) @@ -488,9 +493,9 @@ function _solve!( if iter % nout == 0 && iter > 1 cont += 1 for (norm_Ri, Ri) in zip((norm_Rx, norm_Ry, norm_Rz), @residuals(stokes.R)) - push!(norm_Ri, maximum(abs.(Ri))) + push!(norm_Ri, norm_mpi(Ri) / length(Ri)) end - push!(norm_∇V, maximum(abs.(stokes.R.RP))) + push!(norm_∇V, norm_mpi(stokes.R.RP) / length(stokes.R.RP)) err = max(norm_Rx[cont], norm_Ry[cont], norm_Rz[cont], norm_∇V[cont]) push!(err_evo1, err) push!(err_evo2, iter) diff --git a/test/test_boundary_conditions3D.jl b/test/test_boundary_conditions3D.jl new file mode 100644 index 00000000..9f960392 --- /dev/null +++ b/test/test_boundary_conditions3D.jl @@ -0,0 +1,84 @@ +@static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" + using AMDGPU +elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" + using CUDA +end + +using JustRelax, JustRelax.JustRelax3D +using Test + +const backend = @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU" + AMDGPUBackend +elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA" + CUDABackend +else + CPUBackend +end + +@testset "Boundary Conditions" begin + if backend === CPUBackend + + # test incompatible boundary conditions + @test_throws ErrorException FlowBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=true, front=true, back=true, top=true, bot=true), + ) + + # test with StokesArrays + ni = 5, 5, 5 + stokes = StokesArrays(backend, ni) + stokes.V.Vx .= PTArray(backend)(rand(size(stokes.V.Vx)...)) + stokes.V.Vy .= PTArray(backend)(rand(size(stokes.V.Vy)...)) + stokes.V.Vz .= PTArray(backend)(rand(size(stokes.V.Vz)...)) + + # free-slip + flow_bcs = FlowBoundaryConditions(; + no_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + free_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + ) + flow_bcs!(stokes, flow_bcs) + flow_bcs!(stokes, flow_bcs) # just a trick to pass the CI + + @test @views stokes.V.Vx[ :, :, 1] == stokes.V.Vx[:, :, 2] + @test @views stokes.V.Vx[ :, :, end] == stokes.V.Vx[:, :, end - 1] + @test @views stokes.V.Vx[ :, 1, :] == stokes.V.Vx[:, 2, :] + @test @views stokes.V.Vx[ :, end, :] == stokes.V.Vx[:, end - 1, :] + @test @views stokes.V.Vy[ :, :, 1] == stokes.V.Vy[:, :, 2] + @test @views stokes.V.Vy[ :, :, end] == stokes.V.Vy[:, :, end - 1] + @test @views stokes.V.Vy[ 1, :, :] == stokes.V.Vy[2, :, :] + @test @views stokes.V.Vy[end, :, :] == stokes.V.Vy[end - 1, :, :] + @test @views stokes.V.Vz[ 1, :, :] == stokes.V.Vz[2, :, :] + @test @views stokes.V.Vz[end, :, :] == stokes.V.Vz[end - 1, :, :] + @test @views stokes.V.Vz[ :, 1, :] == stokes.V.Vz[:, 2, :] + @test @views stokes.V.Vz[ :, end, :] == stokes.V.Vz[:, end - 1, :] + + # no-slip + flow_bcs = FlowBoundaryConditions(; + no_slip = (left=true, right=true, front=true, back=true, top=true, bot=true), + free_slip = (left=false, right=false, front=false, back=false, top=false, bot=false), + ) + flow_bcs!(stokes, flow_bcs) + + (; Vx, Vy, Vz) = stokes.V + @test sum(!iszero(Vx[1 , i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 + @test sum(!iszero(Vx[end, i, j]) for i in axes(Vx,2), j in axes(Vx,3)) == 0 + @test sum(!iszero(Vy[i, 1, j]) for i in axes(Vy,1), j in axes(Vy,2)) == 0 + @test sum(!iszero(Vy[i, end, j]) for i in axes(Vy,1), j in axes(Vy,2)) == 0 + @test sum(!iszero(Vz[i, j, 1]) for i in axes(Vz,1), j in axes(Vz,3)) == 0 + @test sum(!iszero(Vz[i, j, end]) for i in axes(Vz,1), j in axes(Vz,3)) == 0 + @test @views Vx[ :, 1, :] == -Vx[ :, 2, :] + @test @views Vx[ :, end, :] == -Vx[ :, end - 1, :] + @test @views Vx[ :, :, 1] == -Vx[ :, :, 2] + @test @views Vx[ :, :, end] == -Vx[ :, :, end - 1] + @test @views Vy[ 1, :, :] == -Vy[ 2, :, :] + @test @views Vy[end, :, :] == -Vy[end - 1, :, :] + @test @views Vy[ :, :, 1] == -Vy[ :, :, 2] + @test @views Vy[ :, :, end] == -Vy[ :, :, end - 1] + @test @views Vz[ :, 1, :] == -Vz[ :, 2, :] + @test @views Vz[ :, end, :] == -Vz[ :, end - 1, :] + @test @views Vz[ 1, :, :] == -Vz[ 2, :, :] + @test @views Vz[end, :, :] == -Vz[end - 1, :, :] + else + @test true === true + end +end \ No newline at end of file