Handling boundaries

[omlins]

Handling boundaries

[ChrisRackauckas]

Okay cool, thought that might be the case.

[ChrisRackauckas]

Ahh yes parallel_indices is the kind of primitive I was looking for, thanks.

[luraess]

Thanks @ChrisRackauckas for reporting and happy to see you could combine ParallelStencil to OrdinaryDiffEq!
Regarding the BCs, there is indeed no specific API available in ParallelStencil. The current example assumes implicitly that T on the boundaries is kept to its initial value, i.e. 0. Imposing Dirichlet BCs in x-direction and Neumann (no flux) in y-direction, leaving the z-direction to 0, one could write

@parallel diffusion3D_step!(du, u, Ci, lam, dx, dy, dz)
T[1  ,:  ,:  ] .= 4 # Dirichlet with value = 4
T[end,:  ,:  ] .= 3 # Dirichlet with value = 3
T[:  ,end,:  ] .= T[:,end-1,:] # Dirichlet
T[:  ,1  ,:  ] .= T[:,2,:] # Neumann no flux

or combine them in a more compact way (and eventually put them into kernels - currently needed for to enable @hide_communication to work with multi-GPUs).

-- EDIT --
When performance is relevant, BCs should be wrapped into a @parallel_indices call:

@parallel_indices (iy,iz) function bc_x!(A::Data.Array)
    A[1  , iy, iz] = 4
    A[end, iy, iz] = 3
    return
end

@parallel_indices (ix,iz) function bc_y!(A::Data.Array)
    A[ix, 1  , iz] = A[ix, 2    , iz]
    A[ix, end, iz] = A[ix, end-1, iz]
    return
end

and called with the according range:

@parallel diffusion3D_step!(du, u, Ci, lam, dx, dy, dz)
@parallel (1:size(du,2), 1:size(du,3)) bc_x!(du)
@parallel (1:size(du,1), 1:size(du,3)) bc_y!(du)

as in e.g. this 3D example.

[luraess]

My example was 2D while yours 3D - see EDITs, also for the kernel version.