Apologies if this topic has been covered or if there is a good guide for this (I have been looking for a while, including at source code for more advanced packages but I can’t wrap my head around if/how they do this sort of thing). I’m not “new” to using Julia (been using it for a few years), but I am new to using Julia “properly”, i.e. writing good code instead of just solving small problems.
I am writing a finite difference code for 2D problems, ideally that other people would be able to use. At the moment I have a loop that computes the second order FD stencil in 2D by looping over the internal nodes:
function SecondDerivative!(uₓₓ::Array,u::Array,Δx::Real,Δy::Real,Index::CartesianIndices)
offx = CartesianIndex(1,0)
offy = CartesianIndex(0,1)
for I in Index
uₓₓ[I] =
(u[I-offx] + 2*u[I] + u[I+offx])/(2.0Δx^2) +
(u[I-offy] + 2*u[I] + u[I+offy])/(2.0Δy^2)
end
uₓₓ
end
The way I’ve written it so far is that when the solver is called, I generate a new function,
D_xx(cache,U) = SecondDerivative!(cache,U,dx,dy,Index)
for calling inside the time loop (since the grid size and number of nodes doesn’t change). But I won’t always want to compute a 5 point stencil. What if I want to be able to use a 2nd order in x and a 4th order in y? For instance I might have these two functions that I want to combine:
function SecondDerivative_Order2(u,Δx)
return (u[j-1] - 2u[j] + u[j+1])/Δx^2
end
function SecondDerivative_Order4(u,Δx)
(1/12*u[j-2] + -4/3*u[j-1] + 5/2*u[j] - 4/3u[j+1] + 1/12*u[j+2])/Δx^2
end
This seems like a job for metaprogramming but I can’t work out how this actually works to combine code snippets.