Skip to content

Commit

Permalink
Add missing 3D no-slip boundary condition (#187)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
albert-de-montserrat authored Jul 4, 2024
1 parent 9e6e0df commit 0682fc9
Show file tree
Hide file tree
Showing 7 changed files with 208 additions and 38 deletions.
1 change: 1 addition & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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[
Expand Down
36 changes: 36 additions & 0 deletions docs/src/man/boundary_conditions.md
Original file line number Diff line number Diff line change
@@ -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
)
```
24 changes: 17 additions & 7 deletions src/Interpolations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
7 changes: 6 additions & 1 deletion src/boundaryconditions/BoundaryConditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))

Expand Down
77 changes: 53 additions & 24 deletions src/boundaryconditions/no_slip.jl
Original file line number Diff line number Diff line change
@@ -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
17 changes: 11 additions & 6 deletions src/stokes/Stokes3D.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
84 changes: 84 additions & 0 deletions test/test_boundary_conditions3D.jl
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 0682fc9

Please sign in to comment.