Bugs in swarm.advection with order=2

[NengLu]

Describe the bug
The particles did not move as expected when using the swam.advection function with order=2.
(run in serial using a manually specified velocity field)

What version of underworld3 are you running and how to reproduce the bug
version of underworld3: lastest dev (in python3.12.7 with petsc3.22.2)
Can reproduce this issue with:

import petsc4py
from petsc4py import PETSc
import underworld3 as uw
from underworld3 import function
import numpy as np
import sympy

dy = 1/4
amplitude = 0.25
offset = 0.     
wavelength = 2
k = 2.0 * np.pi / wavelength

#mdegree,mcont,mradius  = 2,False,dy
mdegree,mcont,mradius = 1,True,dy

mesh = uw.meshing.StructuredQuadBox(elementRes=(8, 8), minCoords=(-1, -1), maxCoords=(1, 1))    
#mesh = uw.meshing.UnstructuredSimplexBox(cellSize=dx,  minCoords=(xmin, ymin), maxCoords=(xmax, ymax),regular=False,refinement=0)

v = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=1,continuous=True)
p = uw.discretisation.MeshVariable("P", mesh, 1, degree=0,continuous=False)

swarm = uw.swarm.Swarm(mesh)
material = uw.swarm.IndexSwarmVariable("M", swarm, indices=2,radius=mradius, proxy_degree=mdegree,proxy_continuous=mcont,update_type=0,npoints=4)  
swarm.populate_petsc(fill_param=3,layout=uw.swarm.SwarmPICLayout.GAUSS)

lightIndex = 0
denseIndex = 1
with swarm.access(material):
    perturbation = offset #+ amplitude * np.cos(k * swarm.particle_coordinates.data[:, 0])
    material.data[:, 0] = np.where(swarm.particle_coordinates.data[:, 1] < perturbation, lightIndex, denseIndex)

x_sym,y_sym = mesh.X[0],mesh.X[1]
vyfn = amplitude * sympy.cos(k * x_sym)*(y_sym+1)/2

with mesh.access(v):
    v.data[:,0] = 0
    v.data[:,1] = uw.function.evaluate(vyfn,v.coords)    
    vmax = v.data[:, 1].max()

step      = 0
time      = 0
dt        = 0
max_time  = 0.5
dt_set    = 1/4

while time < max_time:    
    dt = dt_set
    swarm.advection(V_fn=v.sym, delta_t=dt,order=2,evalf=False)

    courant = vmax * dt / mesh.get_min_radius()
    print(courant)
    
    step += 1
    time += dt

def plot_meshswarm(title,mesh,swarm,material,showSwarm=False,showFig=True,showVel=False,point_size=2,vel=None,skip=2,mag=4,):
    import numpy as np
    import pyvista as pv
    import vtk
    
    pv.global_theme.background = "white"
    pv.global_theme.window_size = [500, 500]
    pv.global_theme.jupyter_backend = "static"
    pv.global_theme.smooth_shading = True
    pv.global_theme.show_edges = True
    pv.global_theme.axes.show = True

    mesh.vtk("tmp_box_mesh.vtk")
    pvmesh = pv.read("tmp_box_mesh.vtk")
    pl = pv.Plotter()
    pl.add_mesh(pvmesh,'Black', 'wireframe')

    if showSwarm:
        with swarm.access():
            points = np.zeros((swarm.data.shape[0],3))
            points[:,0] = swarm.data[:,0]
            points[:,1] = swarm.data[:,1]
        point_cloud = pv.PolyData(points)
        with swarm.access():
            point_cloud.point_data["M"] = material.data.copy()
        pl.add_points(point_cloud,cmap="coolwarm_r",scalars='M',render_points_as_spheres=True, point_size=point_size, opacity=0.8)
    pl.add_title(title,font_size=11)
    pl.remove_scalar_bar()

    if showVel:
        v_vectors = np.zeros((mesh.data[::skip].shape[0], 3))
        v_vectors[:, 0] = uw.function.evaluate(vel.sym[0], mesh.data[::skip], mesh.N)
        v_vectors[:, 1] = uw.function.evaluate(vel.sym[1], mesh.data[::skip], mesh.N)
        v_max = v_vectors.max()
        v_vectors = v_vectors/v_max
        
        v_points = np.zeros((mesh.data[::skip].shape[0], 3))
        v_points[:,0] = mesh.data[::skip][:,0]
        v_points[:,1] = mesh.data[::skip][:,1]
        arrows = pl.add_arrows(v_points,v_vectors, color="black",mag=mag,opacity=1.0, show_scalar_bar=False)
    if showFig:
        pl.show(cpos="xy")
    
    pvmesh.clear_data()
    pvmesh.clear_point_data()

plot_meshswarm('material_advection_order=2',mesh,swarm,material,showSwarm=True,showFig=True,point_size=6)