Leopart Interpolate - cumulative strain tracking

1

Dear all,

For my application I am trying to compute “strain-histories” of particles inside a fluid domain.
I am using LEoPart to initiate and track the particles inside the fluid domain.

Using the particle.interpolate() function, I can compute the instantaneous strain that every particle is exposed to for every timestep.
Now, because I am interested in the accumulation of strain “experience” for the particles, I would like to sum the strain values for every timestep. Using ‘particle.interpolate()’, the new values are used to overwrite the old values, whereas I would like to add the new values to the old values.

The cpp code for interpolate() is as follows:

void particles::interpolate(const Function& phih,
                            const std::size_t property_idx)
{
  std::size_t space_dimension, value_size_loc;
  space_dimension = phih.function_space()->element()->space_dimension();
  value_size_loc = 1;
  for (std::size_t i = 0; i < phih.function_space()->element()->value_rank();
       i++)
    value_size_loc *= phih.function_space()->element()->value_dimension(i);

  if (value_size_loc != _ptemplate[property_idx])
    dolfin_error("particles::interpolate", "get property idx",
                 "Local value size mismatches particle template property");

  for (CellIterator cell(*(_mesh)); !cell.end(); ++cell)
  {
    std::vector<double> coeffs;
    Utils::return_expansion_coeffs(coeffs, *cell, &phih);
    for (std::size_t pidx = 0; pidx < num_cell_particles(cell->index()); pidx++)
    {
      Eigen::MatrixXd basis_mat(value_size_loc, space_dimension);
      Utils::return_basis_matrix(basis_mat.data(), x(cell->index(), pidx),
                                 *cell, phih.function_space()->element());

      Eigen::Map<Eigen::VectorXd> exp_coeffs(coeffs.data(), space_dimension);
      Eigen::VectorXd phi_p = basis_mat * exp_coeffs;

      // Then update
      Point phi_point(_ptemplate[property_idx], phi_p.data());
      _cell2part[cell->index()][pidx][property_idx] = phi_point;
    }
  }
}

I have no experience in cpp, but this might be the moment to start learning.
However, looking at this section of code I believe I need to modify this section:

      // Then update
      Point phi_point(_ptemplate[property_idx], phi_p.data());
      _cell2part[cell->index()][pidx][property_idx] = phi_point;

Instead of “=phi_point”, I would like to do something like “old value”+phi_point.
Using compound assignment, would the following code work according to my needs?
// Then update
Point phi_point(_ptemplate[property_idx], phi_p.data());
_cell2part[cell->index()][pidx][property_idx] += phi_point;

If not, I will address this using python and in the meantime I will be learning cpp.

If yes, how can I compile the c++ code after making the edit?

Thank you in advance,
With kind regards,

2

Does the following link help?

3

Thank you for your quick response and suggestion!

It is certainly helpful. I found this forum post and I read your suggestion in the thread that it could become even more efficient if the code were written in C++.
For this reason, I started exploring the c++ classes and I thought my application might benefit if I could manipulate the “interpolate”-method" in the “particles”-class. Adding the new value to the old value (instead of overwriting the old value), I would be able to track history-dependent properties.

Would you say my suggested change to the method would result in the desired functionality?
Or would this not work/would it require more changes to be made?

4

Adapted from the link I posted. I interpolate a function to one data slot. I then interpolate another function to another data slot. I then add those two particle datasets together, replacing the second data slot.

There may be some formatting issues (which I’m sure you can fix) with my code since I cannot work on a local machine currently.

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()

plot_particles(p)

num_properties = p.num_properties()
for c in cells(mesh):
    for pi in range(p.num_cell_particles(c.index())):
        particle_props = list(p.property(c.index(), pi, prop_num)
                              for prop_num in range(num_properties))

        # Leopart stores particle data in 3D dolfin points                                          
        ap = particle_props[a_idx][0]  # a data                                                     
        bp = particle_props[b_idx][0]  # b data                                                     

        # Reset the particle property data with new values                                          
        p.set_property(c.index(), pi, b_idx, Point(ap + bp))

plot_particles(p)

image
image1128×612 279 KB

image
image1132×598 292 KB

5

No formatting issues, this works flawlessly.

In this way, I am able to compute and track the cumulative strain.
In my adaptation of the above, variable “ap” would contain the strain history and “bp” would contain the instantaneous strain at the current timestep. I multiply “bp” with the timestep “dt” and add it to “ap”

p.set_property(c.index(), pi, b_idx, Point(ap + dt* bp))

With a 120k element mesh and running in parallel on 4 cores, the walltime for every timestep looks like this:
Time solving Navier Stokes equations is 14.357073307037354 s
Time spent computing particle advection is 0.03606915473937988 s
Time spent interpolating strain field to particle is 0.014954328536987305 s
Time spent computing cumulative strain is 0.1708965301513672

Based on this, currently I do not see a necessity to implement tracking of history-dependent properties in c++.

Thank you for the solution you provided. I am able to make a lot of progress with your help.

1 Like