How do I make this function run faster?

What this function does:

I have variables defined on a mesh. Given one of the variables (called 'fieldvar'), this function iterates over every cell of the mesh, iterates over every face of a cell, finds the neighbouring cell. Fetches data of that cell and interpolates to a required degree (given by 'order_new'). Finally, assigns one its faces to 'neigb' and returns it.

Problem:

This funtion is called many times and is unusably slow.

Function definition:

function neighb = neighbour2D(fieldvar)

% function neighb = neighbour2D(fieldvar)
% Purpose: Compute the neighbouring values of fieldvar for every cell commensurate to the order
% demanded. Return array equivalent to e.g. Ez(vmapP), i.e. values of the
% neighbouring cells as a flat array, where every cell contributes
% Nfp*Nfaces no. of entries. Here, the Nfp is acc. to order_new.

Globals2D;

neighb=zeros(length(vmapM),1);
count=0;
for k=1:K
    kself1=fI(k,1); kself2=fI(k,2);
    one=ones(1,Nfp(order_new(k)));
    fmask_local = Fmask{order_new(k)};
    for f=1:Nfaces
        % find neighbor
        kneigh = EToE(k,f); fneigh = EToF(k,f);
        knei1 = fI(kneigh,1); knei2 = fI(kneigh,2);
        val=interp{order_new(kneigh),order_new(k)}*fieldvar(knei1:knei2);
        val=val(fmask_local); val=val(:,fneigh);
    
        % reference length of edge
        v1 = EToV(k,f); v2 = EToV(k, 1+mod(f,Nfaces));
        refd = sqrt( (VX(v1)-VX(v2))^2 + (VY(v1)-VY(v2))^2 );
    
        % find volume node numbers of left and right nodes 
        xself = x(kself1:kself2); yself = y(kself1:kself2);
        xneig = x(knei1:knei2); yneig = y(knei1:knei2);
        x1=xself(fmask_local); x1=x1(:,f);
        x2=interp{order_new(kneigh),order_new(k)}*xneig;
        x2=x2(fmask_local); x2=x2(:,fneigh);
        y1=yself(fmask_local); y1=y1(:,f);
        y2=interp{order_new(kneigh),order_new(k)}*yneig;
        y2=y2(fmask_local); y2=y2(:,fneigh);
        x1=x1*one; x2=x2*one; y1=y1*one; y2=y2*one;
    
        % Compute distance matrix
        D = (x1-x2').^2 + (y1-y2').^2;
        
        [row,~]=find(sqrt(abs(D))<NODETOL*refd);
        neighb(count+1:count+length(val(row)))=val(row);
        count=count+length(val(row));
    end
end

return

Here's what the profiler says:

Lines where the most time was spent

Line Number Code Calls Total Time % Time Time Plot

43 neighb(count+1:count+length(va... 776520 2.198 s 15.3%

42 [row,~]=find(sqrt(abs(D))<N... 776520 1.601 s 11.1%

21 val=interp{order_new(kneigh),o... 776520 1.020 s 7.1%

22 val=val(fmask_local); val=val(... 776520 0.982 s 6.8%

36 y2=y2(fmask_local); y2=y2(:,fn... 776520 0.974 s 6.8%

All other lines 7.611 s 52.9%

[darova]

Can you show the mesh? Why do you need to find neighbours for every cell?

[Rik]

Since it is spending a lot of time on line 43, did you make sure neighb is actually large enough?

Also, if you square NODETOL outside of you loop and remove the sqrt when computing refd, you can skip the sqrt in line 42.

[Apurva]

@darova, the mesh is a square with a circular hole at the center. Here's a zoomed in view:

Here, every triangle is a cell that contains variable no. of nodes (the red/black dots). Computations over any given cell require corresponding data of the neighbouring face (to compute flux; I'm solving a system of hyperbolic PDEs). As the no. of nodes in the neighbouring cell can be dissimilar, neighbouring face data is interpolated accordingly.

[Apurva]

@Rik, I've managed to avoid the sqrt computation as you suggested and that helped. 20 runs now take about 1.5 s lesser. Thank you!

neigb is actually one of the largest vars I have. Every cell has 3 neighbours (except boundaries) and every neighbour contributes some data, so yes, neigb is large. I don't think I fully understand your question about it being actually large enough. How large is justifiable?

[Rik]

The point is that you are storing data in that vector. If it is not large enough yet, Matlab will have to copy the array to extend it, which is very slow.

A=zeros(1,3);
A(1:4)=1:4;%A was too small

[Apurva]

Oh, got it. Yes, I've prealloacted it to the correct size. No implicit expansion should be required.

neighb=zeros(length(vmapM),1);