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
