Commit 047e5297 authored by Jannis Teunissen's avatar Jannis Teunissen

Fix sign in stencil for Neumann boundary condition

This matters for the coarse-grid solver when a non-zero Neumann boundary
condition is used. The sign depends on the direction (high/low) of the boundary.
parent 07376790
......@@ -142,7 +142,7 @@ contains
bc_type = af_bc_dirichlet
bc_val = 0.0_dp
case (af_neighb_highx) ! Higher-x direction
bc_type = af_bc_dirichlet
bc_type = af_bc_neumann
bc_val = 1.0_dp
case default
bc_type = af_bc_neumann
......
......@@ -1173,23 +1173,24 @@ contains
stencil(1, lo(1):hi(1)) = &
stencil(1, lo(1):hi(1)) + &
stencil(nb+1, lo(1):hi(1))
bc_to_rhs(:, nb) = stencil(nb+1, lo(1):hi(1)) * box%dr(nb_dim)
bc_to_rhs(:, nb) = -stencil(nb+1, lo(1):hi(1)) * &
box%dr(nb_dim) * af_neighb_high_pm(nb)
stencil(nb+1, lo(1):hi(1)) = 0.0_dp
#elif NDIM == 2
stencil(1, lo(1):hi(1), lo(2):hi(2)) = &
stencil(1, lo(1):hi(1), lo(2):hi(2)) + &
stencil(nb+1, lo(1):hi(1), lo(2):hi(2))
bc_to_rhs(:, nb) = &
pack(stencil(nb+1, lo(1):hi(1), lo(2):hi(2)) * &
box%dr(nb_dim), .true.)
-pack(stencil(nb+1, lo(1):hi(1), lo(2):hi(2)) * &
box%dr(nb_dim), .true.) * af_neighb_high_pm(nb)
stencil(nb+1, lo(1):hi(1), lo(2):hi(2)) = 0.0_dp
#elif NDIM == 3
stencil(1, lo(1):hi(1), lo(2):hi(2), lo(3):hi(3)) = &
stencil(1, lo(1):hi(1), lo(2):hi(2), lo(3):hi(3)) + &
stencil(nb+1, lo(1):hi(1), lo(2):hi(2), lo(3):hi(3))
bc_to_rhs(:, nb) = &
pack(stencil(nb+1, lo(1):hi(1), lo(2):hi(2), lo(3):hi(3)) * &
box%dr(nb_dim), .true.)
-pack(stencil(nb+1, lo(1):hi(1), lo(2):hi(2), lo(3):hi(3)) * &
box%dr(nb_dim), .true.) * af_neighb_high_pm(nb)
stencil(nb+1, lo(1):hi(1), lo(2):hi(2), lo(3):hi(3)) = 0.0_dp
#endif
case default
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment