Smooth particles hydrodynamic implementation in DealII

[Hassan N]

Dear dealII users

I want to solve a series of conservation law using smooth particles hydrodynamics (SPH). I am looking for a parallel solution.  could you please help me to setup my code correctly.

1- The number of unknowns is 11 (linear momentum 3 components and deformation gradient 9 components). In addition to that I have to save stress tensor, smoothing length, pressure and some other state variables at each particles. So the number of variable that I want to store is huge. Should I store them as 'property' at each particle? 

2- I want to have a particle at each vertex on the triangulation. How to create them?

3- How to create and apply constraints?

4- Finally, (the most important question), in a solution  using the SPH method  computation is done on each particle (Target particle) according to data belong to the particles are around the target particle (Neighbor particles) . How to access to the data that belong to the neighbor particle. For simplicity at this moment we can consider that the neighbor particles are just the particles connected to the target particle on the adjacent cell (one layer). 

Beta regards,

Hassan

[Wolfgang Bangerth]

Hassan,

> I want to solve a series of conservation law using smooth particles
> hydrodynamics (SPH). I am looking for a parallel solution.  could you please
> help me to setup my code correctly.
>
> 1- The number of unknowns is 11 (linear momentum 3 components and deformation
> gradient 9 components). In addition to that I have to save stress tensor,
> smoothing length, pressure and some other state variables at each particles.
> So the number of variable that I want to store is huge. Should I store them as
> 'property' at each particle?

Yes. We have certainly done simulations with many more particle properties.

> 2- I want to have a particle at each vertex on the triangulation. How to
> create them?

Take a look at some of the particle tutorials that show how to create
particles at vertex locations.

> 3- How to create and apply constraints?

You'll have to be more specific. Constraints on what?

> 4- Finally, (the most important question), in a solution  using the SPH
> method  computation is done on each particle (Target particle) according to
> data belong to the particles are around the target particle (Neighbor
> particles) . How to access to the data that belong to the neighbor particle.
> For simplicity at this moment we can consider that the neighbor particles are
> just the particles connected to the target particle on the adjacent cell (one
> layer).

You'll have to have data structures that allow you to identify which particles
are close to which others. The usual one is an rtree, for which deal.II has
wrappers.

The only deal.II-based code I know of that does these sorts of things is Bruno
Blais' 'lethe' code. You should take a look at how this is done inside lethe.

Best
W.

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

[blais...@gmail.com]

If I may suggest, at that point, I would ask you what you are trying to specifically achieve?

If you want to build an SPH code from scratch, this is definitely possible to do using deal.II particle classes (esp. with the exchange and update ghosts functionalities), but this will require a lot of work. In that case, it might also be a good idea to look at an SPH code designed from scratch for this (e.g. DualSPHPhysics is a good one)

In our case, we built our DEM code using deal.II because

A) DEM is simpler than SPH (shorter scale interaction)

B) We wanted to couple our DEM to our CFD and deal.II was the perfect vehicle for that.

It is definitely doable to build from scratch, but it will be a massive endeavour to make it fast.

[Hassan N]

Dear Wolfgang and Bruno 

thank you for your answers. I will study the 'lethe' code. 

I want to apply constraints on linear momentum for a set of particles. 

Actually, I want to build an SPH code from scratch. We are working with a total Lagrangian formulation, so position of particles is constant. Additionally.  we are working with just a layer of neighbor particles in each direction. I did't understand how to find neighbor particles, I was looking for a way to access neighbor particles, something like the one that is used for the cell in step-68 "particle->get_surrounding_cell()".

Thanks

Hassan

[blais...@gmail.com]

Dear Hassan,

There is no way to do that directly in Deal.II

Lethe adds two classes to achieve this using a coarse and then a fine contact detection.

https://github.com/lethe-cfd/lethe/blob/master/include/dem/particle_particle_broad_search.h

https://github.com/lethe-cfd/lethe/blob/master/include/dem/particle_particle_fine_search.h

If you are interested, we could eventually try to port them to deal.II, but these are very lethe specific at the moments.

Otherwise, you can use these ideas to generate your contact list. These will work if you use the background mesh as your contact detection grid (which is what we do, it's quite efficient and robust).

[Hassan N]

Dear Bruno

That is what I am looking for. Surely its good to have them in deal.II. Otherwise, could you please describe shortly the algorithm that you used to do that.  

Thanks, 

Hassan

[blais...@gmail.com]

Dear Hassan,

Everything is described there:

https://arxiv.org/abs/2106.09576

The algorithm is pretty simple. For every cell, we find the list of possible neighbours. We define this list by looking at every cell which shares a vertex with a cell we are presently working with. For all particles, in a cell, we find all potential neighbours using this list of cells which are neighbours of the cell in which the particle lies. That's the coarse contact detection.

The finer one uses a distance function and the radius of the particles to identify if two particles are within potential contact distance (say 1.5 * the particle diameter). Then, these "fine" neighbours are stored in another container and are kept for as many time steps as possible (until we have moved the particle 0.25*diameter). This last part saves a lot of computational time.

The tricky part was choosing the appropriate container for these neighbours because in DEM you need to store the history of a contact. Consequently, you need a data structure that you can easily insert in and remove out of. That's why we went with an unordered map approach. It is maybe not optimal, but it really does work well right now.