Help with ParticleHandler

[Andrew Davis]

I have a question about how to set the quantities for dealii::Particles that are stored in a dealii::ParticleHandler. I have successfully created a dealii::ParticleHandler and randomly placed particles. I have also successfully solved a PDE, whose solution is stored in the vector solution. Then I (again, successfully) interpolate the solution from the PDE onto the particle locations.

Here is the code:

// create the particle handler

dealii::Particles::ParticleHandler<dim> particleHandler(triangulation, mapping, ncomponents);

// randomly place particles in the domain based on some defined pdf 

dealii::Particles::Generators::probabilistic_locations(triangulation, pdf, true, nparticles, particleHandler mapping);

// map the conserved properties onto the particles

dealii::Vector<double> interpolatedParticleQuantities(ncomponents*nparticles);

dealii::Particles::Utilities::interpolate_field_on_particles(dofHandler, particleHandler, solution, interpolatedParticleQuantities);

The problem that I have now is: how do I attach the correct components of interpolatedParticleQuantities with the property of the corresponding particle?

I have tested the implementation by checking:

unsigned int part = 0;

for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); ++iter, ++part ) {

  std::cout << "particle " << part << " quantities: ";

  for( unsigned int i=0; i<ncomponents; ++i ) {

      std::cout << interpolatedParticleQuantities[part*ncomponents + i] << " ";

  }

  std::cout << std::endl;

}

and I have gotten what I expect.  I have also tried attaching the particle properties using: 

unsigned int npart = 0;

for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); ++iter, ++part ) {

  dealii::ArrayView<double> data(interpolatedParticleQuantities.begin() + part*ncomponents, ncomponents);

  iter->set_properties(data);

}

but I get the compiler error:

error: cannot convert ‘dealii::ArrayView<double>’ to ‘const std::vector<double, std::allocator<double> >&’

  353 |       iter->set_properties(data);

Does anyone know what I am doing wrong? Or is there is a better way to do this? I'm sure this is a dumb question, but I couldn't find anything in the tutorials/examples.

[Andrew Davis]

Sorry for the second post---it might also be worth noting that 

unsigned int part = 0;

for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); ++iter, ++part ) {

  std::vector<double> quantities(ncomponents);

  for( unsigned int i=0; i<ncomponents; ++i ) {

    quantities[i] = interpolatedParticleQuantities[part*ncomponents + i];

  }

  iter->set_properties(quantities);

}

results in this exceptions:

The violated condition was: 

    property_pool != nullptr

Additional information: 

    This exception -- which is used in many places in the library -- usually indicates that some condition which the author of the code thought must be satisfied at a certain point in an algorithm, is not fulfilled. An example would be that the first part of an algorithm sorts elements of an array in ascending order, and a second part of the algorithm later encounters an element that is not larger than the previous one.

There is usually not very much you can do if you encounter such an exception since it indicates an error in deal.II, not in your own program. Try to come up with the smallest possible program that still demonstrates the error and contact the deal.II mailing lists with it to obtain help.

[Wolfgang Bangerth]

On 6/16/20 1:46 PM, Andrew Davis wrote:
>
> and I have gotten what I expect.  I have also tried attaching the particle
> properties using:
>
> unsigned int npart = 0;
> for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); ++iter,
> ++part ) {
>   dealii::ArrayView<double> data(interpolatedParticleQuantities.begin() +
> part*ncomponents, ncomponents);
>   iter->set_properties(data);
> }
>
> but I get the compiler error:
>

> *error: *cannot convert ‘*dealii::ArrayView<double>*’ to ‘*const
> std::vector<double, std::allocator<double> >&*’
>
> 353 | iter->set_properties(*data*);

>
>
> Does anyone know what I am doing wrong? Or is there is a better way to do
> this? I'm sure this is a dumb question, but I couldn't find anything in the
> tutorials/examples.

You've already solved this one problem by converting the ArrayView to a
std::vector by hand in your next mail.

> The violated condition was:
>
> property_pool != nullptr
>
> Additional information:
>
> This exception -- which is used in many places in the library -- usually

That's more complicated. The exception happens in the line
iter->set_properties(...)
and tells you that the particle 'iter' points to isn't associated with a
PropertyPool, i.e., the object that stores particle properties. I'm not sure
why that is so. I believe that from your code, you set up the ParticleHandler
to hold properties, but I can't see where the ParticleHandler's PropertyPool
would be propagated to the newly added particles.

Can you create a minimal program that shows the issue? It doesn't have to do
anything useful, just illustrate the error.

Best
W.




--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Andrew Davis]

Thanks for you reply. Here is a minimal example that does 3 three things: (i) interpolate a solution (some simple function in this case, but would actually be the result of solving a PDE) onto the DOFs of a mesh, (ii) interpolate the solution onto randomly placed particles, and (iii) try to associated the interpolated solution with particle properties.  For me, this code compiles but does not run. Do you know why?

I have two other related questions.

(i) In the code below I could replace the counter part with iter->get_id(). Is this the intended use of the ID or would we expect this to fail in more complicated examples where particles may be added/removed or distributed across processors?

(ii) Is it necessary to copy the particle properties into an std::vector or can we just give each particle an iterator the the interolatedParticleQuantities vector? The copy step seems unnecessary and potentially slow/burdensome if there are lots of particles.

#include <deal.II/grid/grid_generator.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/mapping_q1.h>

#include <deal.II/fe/fe_system.h>

#include <deal.II/fe/fe_q.h>

#include <deal.II/numerics/vector_tools.h>

#include <deal.II/particles/generators.h>

#include <deal.II/particles/particle_handler.h>

#include <deal.II/particles/utilities.h>

using namespace dealii;

/// The "solution" to some PDE

template<unsigned int dim>

class SolutionFunction : public Function<dim> {

public:

  inline SolutionFunction() : Function<dim>(dim+1) {}

  inline virtual void vector_value(Point<dim> const& p, Vector<double>& values) const override {

    for( unsigned int i=0; i<dim; ++i ) {

      values[i] = p[i];

    }

    values[dim] = 1.0;

  }

private:

};

int main() {

  // spatial dimension

  constexpr unsigned int dim = 2;

  // number of "unknowns"

  const unsigned int ncomponents = dim+1;

  // number of particles to randomly place

  const unsigned int nparticles = 10;

  // create the mesh

  const MappingQ1<dim> mapping;

  Triangulation<dim> triangulation;

  DoFHandler<dim> dofHandler(triangulation);

  GridGenerator::hyper_cube(triangulation, 0.0, 1.0);

  triangulation.refine_global(5);

  const dealii::FESystem<dim> fe(dealii::FE_Q<dim>(1), ncomponents);

  dofHandler.distribute_dofs(fe);

  // "solve" the pde

  SolutionFunction<dim> function;

  Vector<double> solution(dofHandler.n_dofs());

  VectorTools::interpolate(dofHandler, function, solution);

  // create the particle handler

  Particles::ParticleHandler<dim> particleHandler(triangulation, mapping, ncomponents);

  // randomly place particles in the domain

  ConstantFunction<dim> pdf(1.0);

  Particles::Generators::probabilistic_locations(triangulation, pdf, true, nparticles, particleHandler, mapping);

  // map the conserved properties onto the particles

  dealii::Vector<double> interpolatedParticleQuantities(ncomponents*nparticles);

  Particles::Utilities::interpolate_field_on_particles(dofHandler, particleHandler, solution, interpolatedParticleQuantities);

  // check to make sure we get what we expect

  unsigned int part = 0;

  for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); ++iter, ++part ) {

    std::cout << "particle " << part << " quantities: ";

    for( unsigned int i=0; i<ncomponents; ++i ) {

      std::cout << interpolatedParticleQuantities[part*ncomponents + i] << " ";

    }

    std::cout << std::endl;

    std::cout << "expected quantities: ";

    for( unsigned int i=0; i<dim; ++i ) {

      std::cout << iter->get_location()[i] << " ";

    }

    std::cout << " 1" << std::endl;

    std::cout << std::endl;

  }

  // now try to set the interpolated quantities as particle properties

  part = 0;

  for( auto iter=particleHandler.begin(); iter!=particleHandler.end(); ++iter, ++part ) {

    std::vector<double> quantities(ncomponents);

    for( unsigned int i=0; i<ncomponents; ++i ) {

      quantities[i] = interpolatedParticleQuantities[part*ncomponents + i];

    }

    iter->set_properties(quantities);

  }

}

[Wolfgang Bangerth]

On 6/18/20 9:42 AM, Andrew Davis wrote:
> Thanks for you reply. Here is a minimal example that does 3 three things: (i)
> interpolate a solution (some simple function in this case, but would actually
> be the result of solving a PDE) onto the DOFs of a mesh, (ii) interpolate the
> solution onto randomly placed particles, and (iii) try to associated the
> interpolated solution with particle properties.  For me, this code compiles
> but does not run. Do you know why?

Yes, this is a bug. I've opened
https://github.com/dealii/dealii/issues/10584
Let's take the discussion of the failure there. I think I've got a fix.

> (i) In the code below I could replace the counter part with iter->get_id(). Is
> this the intended use of the ID or would we expect this to fail in more
> complicated examples where particles may be added/removed or distributed
> across processors?

It depends on what you want to do. Your counter goes over the locally owned
particles and takes their index *within the locally owned particles*.
iter->get_id() returns an id that is intended to be globally unique across all
processors.

> (ii) Is it necessary to copy the particle properties into an std::vector or
> can we just give each particle an iterator the the
> interolatedParticleQuantities vector? The copy step seems unnecessary and
> potentially slow/burdensome if there are lots of particles.

You mean in this part of the code?

for( auto iter=particleHandler.begin(); iter!=particleHandler.end();
++iter, ++part ) {
std::vector<double> quantities(ncomponents);
for( unsigned int i=0; i<ncomponents; ++i ) {
quantities[i] = interpolatedParticleQuantities[part*ncomponents + i];
}

iter->set_properties(quantities);
}

Yes, this was an oversight that I fixed a while ago already (but post 9.2):
ParticleAccessor::set_properties() did not have an overload for ArrayView.
With the current development version, you can write this as follows:

for( auto iter=particleHandler.begin(); iter!=particleHandler.end();
++iter, ++part )

iter->set_properties
(make_array_view(interpolatedParticleQuantities.begin()+
part*ncomponents,
interpolatedParticleQuantities.begin()+
part*ncomponents + ncomponents));