[Leopart] Update particle properties manually

1

Hello
I am currently using Fenics and Leopart to model the evolution of ice crystal orientation in a stokes velocity field. I would like to be able to manually update the properties carried by the Leopart particles at each timestep using a model I have implemented in numpy. At the moment I have this working by extracting the properties, modifying them and then creating new particles (with the same location) at each timestep. However I was wondering if it is possible to more elegantly assign new values to the particle properties. This approach would I assume be more suited to parallelization.

I’ve attached a minimal example of how I’m currently updating particle properties:

from fenics import*
from mshr import *
import leopart as lp
import numpy as np


nt=10
dt = 0.01


npart = 15
npmin=12
npmax=20


mesh = UnitSquareMesh(20,20)

#Velocity 
V = VectorFunctionSpace(mesh, "DG", 1)
P = FunctionSpace(mesh,"DG",1)
u = Function(V)
u.assign(Expression(("x[1]","0.0"),degree=1))
(ux,uy)=split(u)



# Initialize particles
x = lp.RandomCell(mesh).generate(npart)

def numpy2particle(xnew,q,r):
    return lp.particles(xnew,[q, r], mesh)


def particle2numpy(p):
    x = p.positions()
    a = np.array(p.get_property(1))
    b = np.array(p.get_property(2))
    return x,a,b

# Initial particle
a = np.zeros(x.shape[0])
b = np.zeros(x.shape[0])
p = numpy2particle(x,a,b)

# Functions for interpolating
a_df = Function(P)
b_df = Function(P)

for it in range(nt):

     #Advect particles
    ap = lp.advect_rk3(p, V, u, "open")
    ap.do_step(dt)

    #Add delete
    AD = lp.AddDelete(p,npmin,npmax,[a_df,b_df])
    AD.do_sweep()


    x,a,b = particle2numpy(p)

    ## Change properties
    a= 2*b
    b= 5. + a        ##just as a minimal example 

    # Update particle by creating a new particle
    p = numpy2particle(x,a,b)

    # Create fenics functions for add delete
    lstsqa = lp.l2projection(p,P,1)
    lstsqa.project(a_df)

    lstsqb = lp.l2projection(p,P,1)
    lstsqb.project(b_df)

Thanks for the help

1 Like
2

Consider the following which directly manipulates each particle belonging to a provided list of candidate cells:

from dolfin import *
import leopart as lp
import numpy as np

npart = 15

mesh = UnitSquareMesh(20, 20)
P = FunctionSpace(mesh, "DG", 1)

# Initialize particles
x = lp.RandomCell(mesh).generate(npart)

# Initial particle
a = np.zeros(x.shape[0])
b = np.zeros(x.shape[0])

# Property indices
a_idx, b_idx = 1, 2

p = lp.particles(x, [a, b], mesh)

# Some initial data
a_df_0 = interpolate(Expression("x[0]", degree=1), P)
b_df_0 = interpolate(Expression("x[1]", degree=1), P)
p.interpolate(a_df_0, a_idx)
p.interpolate(b_df_0, b_idx)


def plot_particles(particles):
    positions = particles.positions()
    import matplotlib.pyplot as plt
    idxs_to_plot = (a_idx, b_idx)
    fig, axs = plt.subplots(1, len(idxs_to_plot))
    for prop_idx in idxs_to_plot:
        axs[prop_idx-1].scatter(positions[:,0], positions[:,1],
                                c=particles.get_property(prop_idx), s=0.1)
        axs[prop_idx-1].set_aspect("equal")
        axs[prop_idx-1].set_title(f"particle property {prop_idx}")
    plt.show()


def edit_properties(particles, candidate_cells):
    num_properties = particles.num_properties()
    num_particles_changed = 0

    for c in candidate_cells:
        for pi in range(particles.num_cell_particles(c)):
            particle_props = list(particles.property(c, pi, prop_num)
                                  for prop_num in range(num_properties))

            # Leopart stores particle data in 3D dolfin points
            px = particle_props[0]  # position is always idx 0
            ap = particle_props[a_idx][0]  # a data
            bp = particle_props[b_idx][0]  # b data
            # Property idxs > b_idx are related to the RK advection scheme

            # Reset the particle property data with new values
            particles.set_property(c, pi, a_idx, Point(2 * bp))
            particles.set_property(c, pi, b_idx, Point(5 + ap))

            # Record the number of particles changed
            num_particles_changed += 1

    info(f"Changed {num_particles_changed} particles")


# We change all particles in all cells
all_cells = [c.index() for c in cells(mesh)]

# Direct manipulation method
plot_particles(p)
edit_properties(p, all_cells)
plot_particles(p)

Before editing properties:
image
After editing properties:
image

Bear in mind that if you change the particles’ positions using this direct manipulation method, you should call particles.relocate() to redistribute the data across cells and MPI processes.

It is also likely much more efficient to write the above code in C++. The loops are not very efficient in python. But are “good enough” if you don’t need to edit all the particles in the mesh.

For context, this is the basis for the numerical implementation of the melting model in this work.

1 Like
3

Hi Nate, thanks for your quick and helpful response. Unfortunately I do need to edit all the particles in the mesh. Despite this the python loops are of similar speed to my current implementation, so a c++ implementation would likely be significantly faster than what I have. I don’t really have much experience with c++, what do you think would be the best way to approach this? Would I have to add something into leopart or could I do something with pybind?

I also had to modify your code slightly as my minimal implementation was perhaps too minimal, In my problem I assign the properties from a numpy vector calculated externally rather than just doing ‘b=5+a’