omlins / ParallelStencil.jl

Package for writing high-level code for parallel high-performance stencil computations that can be deployed on both GPUs and CPUs
BSD 3-Clause "New" or "Revised" License
311 stars 31 forks source link

Question about @inn macro #61

Closed ttaczak closed 2 years ago

ttaczak commented 2 years ago

I am relatively new to metaprogramming in Julia so this may simply be a personal issue.

With that acknowledged, I am struggling to understand how the @inn macro in FiniteDifferences.jl behaves as advertised. The documentation defines @doc "@inn(A): Select the inner elements of A. Corresponds to A[2:end-1]" :(@inn). The relevant definitions in the file include const ix = INDICES[1] const ixi = :($ix+1) ... macro all(A::Symbol) esc(:( $A[$ix ] )) end macro inn(A::Symbol) esc(:( $A[$ixi ] )) end From this definition, it seems that @inn would simply shift all of the indices of A::Symbol by one value (resulting in an out of bounds error). Considering that all of the miniapps function as intended, I know I must be missing something but I can't for the life of me figure out how it works.

luraess commented 2 years ago

Hi @ttaczak ,

Thanks for your interest in using ParallelStencil!

Note: Since your question is not an issue with the package but rather an application / usage question, it would be better suited to post it on Julia Discourse in the most appropriate section as suggested here - maybe in Modelling & Simulations topic ? As such, other folks can also profit from the discussion and we keep GitHub issues for actual technical issues with the package itself.

Regarding your question:

I am struggling to understand how the 'https://github.com/INN' macro in FiniteDifferences.jl behaves as advertised

As stated from querying the doc via the help ? within the REPL:

julia> using ParallelStencil

julia> using ParallelStencil.FiniteDifferences

julia> using ParallelStencil.FiniteDifferences2D

help?> @inn
  @inn(A): Select the inner elements of A. Corresponds to A[2:end-1,2:end-1].

help?> @all
  @all(A): Select all elements of A. Corresponds to A[:,:].

julia>

Which means that @inn(A) would be needed to select the inner elements of A, which corresponds to A[2:end-1,2:end-1] in e.g. 2D vectorised form. This would be used when computing e.g. the second derivatives where one obviously looses the boundary values of the array A as in the following explicit diffusion example (skipping the diffusion coefficient - or having it =1):

@inn(A2) = @inn(A) + dt * ( @d2_xi(A)/dx^2 + @d2_yi(A)/dy^2 )

The vectorised analogy would be

A2[2:end-1,2:end-1] .= A[2:end-1,2:end-1] .+ dt .* ( diff(diff(A[:,2:end-1],dims=1),dims=1)./dx^2 + diff(diff(A[2:end-1,:],dims=2),dims=2)./dy^2 )

Hope this helps!

ttaczak commented 2 years ago

Thank you for your quick response! I believe I was perhaps a bit unclear in my question. I understand how the function works since I have access to its definition.

I am primarily curious why it works as intended because I am trying to create my own set of macros similar to those defined in the FiniteDifference.jl file included in the ParallelStencil package. I will use @all and simply xi an example because I believe these will better illustrate my confusion.

Using the vectorized notation @all(A) would be equivalent to A[1:end]. The @all macro is defined as

macro all(A::Symbol) esc(:( $A[$ix ] )) end

From this, I would assume that a macro which references $A[$ix+1] would be equivalent to A[2:end+1], thereby resulting in an out of bounds error. However, my understanding is demonstrably false, since @d (which works without any errors) is defined as

macro d(A::Symbol) esc(:( $A[$ix+1] - $A[$ix] )) end

luraess commented 2 years ago

I get your point now. I guess the confusion comes from the fact that these macro need some bound check to perform correctly. This is performed implicitly in the @parallel kernel, where we assume that the bounds of the operation are the size of the left-hand.side expression, Here a small example:

nx = 7
A  = rand(nx)
B  = rand(nx-2)

and computing

@inn(A) = @all(B)

would not fail as the loop equivalent would be

for ix=1:nx-2
    A[ix+1] = B[ix]
end

where we know the range being nx-2 as we want inner points. Then we need the +1 shift to have the [2:end-1] behaviour. Alternative would be to have a shifted range

for ix=2:nx-1
    A[ix] = B[ix-1]
end

but this is less handy to generalise. Hope this helps!

ttaczak commented 2 years ago

Thank you so much! I really appreciate your help