FluidityProject / fluidity

Fluidity
http://fluidity-project.org
Other
366 stars 115 forks source link

Active particle scheme #284

Closed Cmath2 closed 3 years ago

Cmath2 commented 4 years ago

This pull request aims to implement an active particle scheme within the fluidity framework. The build extends upon the now merged passive particle scheme, and enables particles to interact with, and influence the values of fields during advection. This is done in three steps:

  1. Relevant field values are interpolated onto particles at their current position in space.
  2. Particles then calculate 'attributes', based on their position, the current simulation time, and interpolated field values.
  3. Particle attributes are then interpolated back onto the nodes of a chosen mesh, with interpolated particle values able to generate new fields, or influence old ones.

With particles able to interact with and generate fields, limits must be placed on the number of particles present in set regions of the domain. This build has two primary components, the implementation of an active particle scheme able to interact with the fields of a simulation, and the implementation of spawning and deleting algorithms to control the number of particles present (per control volume). Currently, particles are spawned by cloning all parent particles that exist within a control volume and perturbing their location by a user specified amount. Conversely, particles are currently deleted by implementing a coin flip algorithm, where a coin is flipped for every particle within a control volume, and individual particles are deleted if the result is heads.

The schema has been updated to include new options for particle spawning and deleting (under each particle_group). Here the user can set the minimum and maximum thresholds of particles per control volume, as well as the location of spawned particles, either by setting a radius around the parent particle in which the new particle can spawn, or by allowing particles to spawn randomly within the current control volume prioritizing elements with fewer particles. Spawning caps can also be placed on which particle subgroups are spawned, where the ratio of particles in a specific group to the total particles in the control volume must be equal to or greater than the cap for spawning of that group to occur (e.g., if a cap of 90% is set, 90% or more of the particles in a control volume must belong to one group for spawning of that group to occur). Finally, the user can specify whether spawned particles will copy attribute values from their parents, or have their attribute values calculated as a weighted function of all particles in the control volume.

Options have also been added underneath diagnostic scalar fields, with two new algorithms, 'from_particles' and 'number_of_particles', now available. -The 'from_particles' algorithm implements the ratio method to interpolate the given particle_group and particle_attribute values onto a user specified mesh (Currently only works for p1 meshes). Two different interpolation options are also present here, being a weighted distance interpolation of particle attributes to nodes, or nearest neighbour interpolation. -The 'number_of_particles' algorithm simply counts the number of particles present per control volume, and sets this as the value for the relevant node. This enables users to track regions of a domain that are relatively dense, or sparse in particles.

Additional smaller changes have also been made, which include: -With the recent h5part IO changes, the option to make particles static or write NaN outputs upon leaving a domain have been removed. Instead particles which leave a domain are simply deleted. -Detector%types have been removed. By default all detectors are now static, and all particles are now lagrangian -Detector%name has been removed, as they are not used in the IO, and detectors/particles are uniquely identified by the combination of their id and proc_id -Write_diagnostics has also been updated to include the total number of particle present per particle_group, as well as each particle_subgroup within the particle_group -It is no longer required to define the number of particles for each particle sub_group in the schema file -Particles can now be initialised after the start of the simulation. This option appears in the schema file as initialise_during_simulation under particle_subgroup -Profiler calls have been added to time various aspects of the particle scheme -Several test cases have been added to test various features of the active particle scheme -Several bug fixes to do with particle parallelization -Updates to test case option files to account for changes to the schema

drhodrid commented 3 years ago

@Cmath2 - since @angus-g has looked through the code carefully and @stephankramer will doubtless do so as well, following my review, I looked at the other aspects. In general, I think you've done a great job of ticking all boxes! The only things missing in my opinion are:

  1. Manual - the manual needs to be updated to highlight what can now be done with "active" particles and the various underlying assumptions/methods documented. In addition, .flml options need documentation in the manual.
  2. Given that this is some blooming awesome new Fluidity functionality, I would recommend adding an example highlighting a commonly used case (under examples). I'd suggest the mu10 RT case for this as it's a little bit awesome. See examples/stokes_square_convection for a suggested structure to follow.
  3. At least 3 of your tests are too long for the short or medium test suite. These are particle_stable_layer, particle_unstable_layer and particle_vortex_3d. These should be moved across to longtests post-merge (merely flagging this now so it's recorded). If there are others that I'm missing, please remember to move those across as well.
  4. The comments in the .rnc/rng are fine in general, but a few could do with updating (e.g. text that comes up when you activate "particle_group".
  5. Test case .xmls - I'd recommend adding a little more information on what is being tested and how it is being tested to the test descriptions. This will ensure that years down the line, if specific tests start to fail we have a better idea of where to look for the problem.
  6. Test case names - perhaps be a little more specific for some? E.g. Particle_multimaterial --> particle_rayleigh_taylor_instability (or similar). Others (e.g. particle_spawning) are as clear as they can be.

Aside from that, you've done blooming awesome. Now over to @stephankramer for the more careful checks!

drhodrid commented 3 years ago

Thanks for doing such a thorough job addressing my comments @Cmath2. I guess we now hand over to @stephankramer for the ultimate test!!

stephankramer commented 3 years ago

Thanks, I had indeed waited for @angus-g - who with @drhodrid, seem to have done a pretty thorough job already - let's see if I can still find something to moan about..

Cmath2 commented 3 years ago

Ok @stephankramer I believe the branch is ready for you again. We had a couple of bugs arise regarding attribute_arrays which have been handled, and I've also added another test case to ensure they work.

Cmath2 commented 3 years ago

I believe everything is up to date and ready for the next round @stephankramer

Cmath2 commented 3 years ago

Hey @stephankramer, I believe I've addressed everything except for the random_seed call.

For this the goal is to create a random seed so there is consistency in the random numbers generated for spawning/deleting between runs. I do see the point that calling random_seed every time particle_cv_check is called is redundant though. If we are to keep the seed, would you suggest I instead set it within Fluids.F90 before the first call to particle_cv_check (e.g. on line 342)? Or somewhere else like in initialise_particles with a check that the spawning/deleting option is active?

Regarding your other feedback, I don't believe I've done any thorough profiling of the current spawning/deleting scheme. It probably would be best to leave this one for Angus in a future commit.

stephankramer commented 3 years ago

Yeah there's probably some bike-shedding we could do about the random seed. One of the potential issues is that it affects other places in the code that uses random numbers - and in some cases you want deterministic behaviour and sometimes you don't (turbulence modelling). Probably the general solution is to have a central place to do this, a user choice whether things are deterministic or not, and also a mechanism to checkpoint the seed value. At the moment I think this only affects synthetic bcs with LES turbulence modelling, which nobody uses (in particular not in combination with particles!) - so let's just ignore this for now.

Cmath2 commented 3 years ago

I believe that covers all of the tabs in the schemas @stephankramer, let me know if I missed any. With that wrapped up I think we should be good to merge when the light goes green :)