Extremely slow field interpolation

dolfinx
1

Hi,

I am solving a 2d transient problem where I need to visualize fields into an animation. At the moment I save the solved field at each time step of interest in a list and later in the post-processing I interpolate fields in a numpy array which is plotted frame by frame. My motivation for doing it this is that I am comfortable with customizing matplotlib plots so it helps presenting things in way I need.

However, the field interpolation is turning out to be extremely slow - it takes roughly around 5 seconds to interpolate the two vector components from the solved scalar potential! Once I solve for large enough number of points it would be very impractical to keep up with this.

I wanted to know how others process the time domain results? Is using paraview and/or pyvista following the tutorial approach more efficient? Or it will be similar?

Thanks a lot

2

Are you compiling a UFL expression before interpolation (through a dolfinx.fem.Expression object)? If so did you make sure to compile it only once and not at each time step? Recompiling it multiple times can slow down the program a lot.

3

Thanks for your quick reply! I think I am indeed doing it the way you suggest. To be sure, here is the code:

OmegaSol = dolfinx.fem.FunctionSpace(self.mesh, ("DG", 2))
Ex, Ey = dolfinx.fem.Function(OmegaSol), dolfinx.fem.Function(OmegaSol)
Power = np.zeros((len(saved_steps), pnts.shape[1]))
for step in tqdm(saved_steps):
       Ex_expr = dolfinx.fem.Expression(-self.V_list[step].dx(0), OmegaSol.element.interpolation_points)
       Ey_expr = dolfinx.fem.Expression(-self.V_list[step].dx(1), OmegaSol.element.interpolation_points)
       Ex.interpolate(Ex_expr)
       Ey.interpolate(Ey_expr)
       Ex_sol = model_A.field_interp('Ex', points)
       Ey_sol = model_A.field_interp('Ey', points)
       Power[step,:] = sigma * (Ex_sol **2 + Ey_sol **2) # J.E #

In the above, ‘self.V_list’ is the list containing the solved potential at chosen steps. Please ignore field_interp() function, it is a wrapper to interpolation code I found from one of the tutorials. But for the approach, is this what you were suggesting in your comments? If so, could you point out where it needs to change with regards to my example?

Thanks!

4

In your code you are creating expressions at each loop iteration. Hence you are calling the JIT compiler 2*saved_steps times, which can be quite slow.

I would suggest creating these expressions before the for loop, and then just perform the interpolation inside the loop.

OmegaSol = dolfinx.fem.FunctionSpace(self.mesh, ("DG", 2))
Ex, Ey = dolfinx.fem.Function(OmegaSol), dolfinx.fem.Function(OmegaSol)
Power = np.zeros((len(saved_steps), pnts.shape[1]))

# create a dummy function to build Ex/Ey expressions
v = dolfinx.fem.Function(self.V_list[0].function_space)
Ex_expr = dolfinx.fem.Expression(-v.dx(0), OmegaSol.element.interpolation_points)
Ey_expr = dolfinx.fem.Expression(-v.dx(1), OmegaSol.element.interpolation_points)

for step in tqdm(saved_steps):
       # copy saved data inside interpolation variable
       v.vector[:] = self.V_list[step].vector[:]

       Ex.interpolate(Ex_expr)
       Ey.interpolate(Ey_expr)
       Ex_sol = model_A.field_interp('Ex', points)
       Ey_sol = model_A.field_interp('Ey', points)
       Power[step,:] = sigma * (Ex_sol **2 + Ey_sol **2) # J.E #

With this code the compiler is called only twice.
Can you try if something like this is faster?

2 Likes
5

I just did, unfortunately it does not offer any noticeable improvement in compute time.

6

It is unclear what these functions do, as you have not supplied a minimal working example reproducing your issue.

Also, it would be very interesting to see the breakdown of what takes time in your code, i.e. how much time does Ex.interpolate and Ey.interpolate actually take in your code?

7

Thanks for your reply.

I agree, while I prepare an MWE, here is what goes in the wrapper model_A.field_interp(‘Ex’, points)

    def field_interp(self, field:str, points:np.ndarray)->np.ndarray:
        field_out = FD.interp2(self.mesh, getattr(self, field), points=points)
        return field_out
# in the file referenced as FD, here is the edfinition of interp2

def interp2(mesh, uh, points:typing.List, ndim = -1)->np.ndarray:
    '''Interpolate the given model over specified points. Adapted from
    https://github.com/jorgensd/dolfinx-tutorial/blob/v0.4.0/chapter1/membrane_code.ipynb
    Args:
        mesh: fenics mesh object
        uh: ufl object, either a solution of the problem or another variable defined on mesh
        points: numpy array to interpolate on, size ndim, num_points
    Output:
        uh_values: np.array of shape (points.shape[1],)
    '''
    if np.shape(points)[0] != 2:
        raise Exception("Points array must has a shape of 2 x num_points for 2D surface interpolation")
    else:
        points = np.array(points) # make sure it is an np array
        points = np.vstack((points, np.zeros((1, points.shape[1]))))

    # first task is to identify all the mesh elements inside the
    bb_tree = dolfinx.geometry.BoundingBoxTree(mesh, mesh.topology.dim)

    cells = []
    points_on_proc = []
    # Find cells whose bounding-box collide with the the points
    cell_candidates = dolfinx.geometry.compute_collisions(bb_tree, points.T)
    # Choose one of the cells that contains the point
    colliding_cells = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
    uh_values = np.zeros(points.shape[1])
    uh_values[:] = np.nan
    for i, point in enumerate(points.T):
        if len(colliding_cells.links(i)) > 0:
            # points_on_proc.append(point)
            # cells.append(colliding_cells.links(i)[0])
            cell = colliding_cells.links(i)[0]
            uh_values[i] = uh.eval(point, cell)
    return uh_values

I will also do time profiling and report.

8

If the points and cells are constant in time, calling eval is highly inefficient.
Instead you should pull the physical points back to the reference cell.

The code you produced is not a minimal working example, I cannot copy this code and reproduce any timings etc.

I would strongly suggest making an example on a unit square or cube, which illustrates the operation you want to perform.

9

Thanks a lot for your feedback and hint! I will be back soon with further information on this.

10

For now, could you refer to some examples based on an efficient alternative to using uh.eval the way I did? Indeed, points and cells are constant in time.