Question:

I’m currently implementing an Immersed Boundary Method (IBM) in Julia, and I’ve written a function to compute interpolation weights for IBM forcing. I want my IBM code to be differentiable for which I will be using Enzyme.jl. My enzyme autodiff function is showing that compilation has failed when I am using this function inside the simulation block. How can I make this function differentiable? Here’s the function:

``````function calculate_interpolation_weights_indices(x_grid::Array{Float64}, y_grid::Array{Float64},
lag_x::Array{Float64}, lag_y::Array{Float64}, n_neigh::Int64)
# n_neigh are the number of neighbours that the point will have in one direction
lag_len=length(lag_x); # Total number of lagrange points
x_len=length(x_grid);  # Total number of eulerian x indices
y_len=length(y_grid);  # Ttoal number of eulerian y indices
points_neigh=(2*n_neigh)^2; # Number of neighbours for each point
# Defining the weights matrices
weights_vel=zeros(Float64,lag_len,points_neigh);
weights_force=zeros(Float64,lag_len,points_neigh);
weights_c_f=zeros(Float64,lag_len,points_neigh);
# Defining the indices matrices
i_ind=zeros(Int64,lag_len,points_neigh);
j_ind=zeros(Int64,lag_len,points_neigh);

# It is assumed that the lagrangian points form a closed curve
# Thus P0-P1-P2.....P(N-1)-P0 is the curve in question
# This is done to compute the arc lengths that will be essesntial for weight calculations

for ind=1:lag_len
x_cord=lag_x[ind]; # x coordinate of the lagrangian point
y_cord=lag_y[ind]; # y coordinate of the lagrangian point
# Now let us find the nearest indices of the points
near_i=1;
near_j=1;
for i=1:(x_len-1)
if (x_cord>=x_grid[i])&&(x_cord<x_grid[i+1])
near_i=i+0;
break
end
end

for j=1:(y_len-1)
if (y_cord>=y_grid[j])&&(y_cord<y_grid[j+1])
near_j=j+0;
break
end
end
# Now we need to precompute some quatitites and then
# Calculate the neighbours and the indices
dx=x_grid[near_i+1]-x_grid[near_i];
dy=y_grid[near_j+1]-y_grid[near_j];
dh=sqrt(dx*dy); # Thickness of the immersed boundary
# Now we need to find the length of the immersed boundary element
fwd=ind+1;
back=ind-1;
if ind==1
back=lag_len;
elseif ind==lag_len
fwd=1;
end
df=sqrt((lag_x[fwd]-lag_x[ind])^2+(lag_y[fwd]-lag_y[ind])^2);
db=sqrt((lag_x[back]-lag_x[ind])^2+(lag_y[back]-lag_y[ind])^2);
ds=0.5*(df+db);
# Now we are ready to find out the interpolation weights
curr=1;
for i=(near_i-n_neigh+1):(near_i+n_neigh)
for j=(near_j-n_neigh+1):(near_j+n_neigh)
x_eu=x_grid[i]; # X Coordinate of the euler point
y_eu=y_grid[j]; # Y Coordinate of the euler point
rx=(abs(x_eu-x_cord))/dx;
ry=(abs(y_eu-y_cord))/dy;
delta_x=0.0;
delta_y=0.0;
if rx<=2.0
delta_x=0.25*(1+cos(pi*0.5*rx));
end
if ry<=2.0
delta_y=0.25*(1+cos(pi*0.5*ry));
end
delta=delta_x*delta_y;
i_ind[ind,curr]=i+0;
j_ind[ind,curr]=j+0;
weights_vel[ind,curr]=delta+0.0;
weights_force[ind,curr]=delta*(ds*dh)/(dx*dy);
weights_c_f[ind,curr]=delta*ds*dh;
curr=curr+1;
end
end

end
return points_neigh,weights_vel,weights_force,weights_c_f,i_ind,j_ind
end
``````

Hi!
Can you share the full reproducible example, showing the way you called Enzyme and the resulting error?

diff_test.jl (3.2 KB)
kt_acm.jl (16.8 KB)
The solver code is kt_acm.jl and a simple test is diff_test.jl

Can you also post the full error message? Otherwise it’s unclear what is going wrong

This is a snippet of the code diff_test.jl.
If i let the x_eu and y_eu both be Const. The execution fails but when I change them to Duplicated, the code hangs and I don’t get any output.

Can you post a text of the full stack trace displayed in the REPL, instead of a screenshot from the text file?

``````warning: didn't implement memmove, using memcpy as fallback which can result in errors
warning: didn't implement memmove, using memcpy as fallback which can result in errors
warning: didn't implement memmove, using memcpy as fallback which can result in errors
warning: didn't implement memmove, using memcpy as fallback which can result in errors
Type tree: {[-1]:Pointer, [-1,0]:Pointer, [-1,0,-1]:Float@double, [-1,8]:Integer, [-1,9]:Integer, [-1,10]:Integer, [-1,11]:Integer, [-1,12]:Integer, [-1,13]:Integer, [-1,14]:Integer, [-1,15]:Integer, [-1,16]:Integer, [-1,17]:Integer, [-1,18]:Integer, [-1,19]:Integer, [-1,20]:Integer, [-1,21]:Integer, [-1,22]:Integer, [-1,23]:Integer, [-1,24]:Integer, [-1,25]:Integer, [-1,26]:Integer, [-1,27]:Integer, [-1,28]:Integer, [-1,29]:Integer, [-1,30]:Integer, [-1,31]:Integer, [-1,32]:Integer, [-1,33]:Integer, [-1,34]:Integer, [-1,35]:Integer, [-1,36]:Integer, [-1,37]:Integer, [-1,38]:Integer, [-1,39]:Integer}
You may be using a constant variable as temporary storage for active memory (https://enzyme.mit.edu/julia/stable/#Activity-of-temporary-storage). If not, please open an issue, and either rewrite this variable to not be conditionally active or use Enzyme.API.runtimeActivity!(true) as a workaround for now

Stacktrace:
[1] getindex
@ ./tuple.jl:29
[2] iterate (repeats 2 times)
@ ./tuple.jl:69
[3] vcat
@ ./array.jl:1812

Stacktrace:
[1] throwerr(cstr::Cstring)
@ Enzyme.Compiler ~/.julia/packages/Enzyme/ZH7JA/src/compiler.jl:2914
in expression starting at /home/arijit/Desktop/Julia res/KT_acm/diff_test.jl:89
``````

This is the result I am getting when I keep x_eu,y_eu as Const.