From: Rene Gassmoeller Date: Tue, 8 Aug 2017 20:01:53 +0000 (-0600) Subject: Initial particle implementation X-Git-Tag: v9.0.0-rc1~1013^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c645548d0d66aacc4554c64b54369d025a1a8310;p=dealii.git Initial particle implementation --- diff --git a/include/deal.II/particle/particle.h b/include/deal.II/particle/particle.h new file mode 100644 index 0000000000..ab488769b3 --- /dev/null +++ b/include/deal.II/particle/particle.h @@ -0,0 +1,309 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__particle_particle_h +#define dealii__particle_particle_h + +#include + +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Particle +{ + /** + * A namespace for all type definitions related to particles. + */ + namespace types + { + /* Type definitions */ + +#ifdef DEAL_II_WITH_64BIT_INDICES + /** + * The type used for indices of particles. While in + * sequential computations the 4 billion indices of 32-bit unsigned integers + * is plenty, parallel computations using hundreds of processes can overflow + * this number and we need a bigger index space. We here utilize the same + * build variable that controls the dof indices because the number + * of degrees of freedom and the number of particles are typically on the same + * order of magnitude. + * + * The data type always indicates an unsigned integer type. + */ + typedef unsigned long long int particle_index; + + /** + * An identifier that denotes the MPI type associated with + * types::global_dof_index. + */ +# define ASPECT_PARTICLE_INDEX_MPI_TYPE MPI_UNSIGNED_LONG_LONG +#else + /** + * The type used for indices of particles. While in + * sequential computations the 4 billion indices of 32-bit unsigned integers + * is plenty, parallel computations using hundreds of processes can overflow + * this number and we need a bigger index space. We here utilize the same + * build variable that controls the dof indices because the number + * of degrees of freedom and the number of particles are typically on the same + * order of magnitude. + * + * The data type always indicates an unsigned integer type. + */ + typedef unsigned int particle_index; + + /** + * An identifier that denotes the MPI type associated with + * types::global_dof_index. + */ +# define ASPECT_PARTICLE_INDEX_MPI_TYPE MPI_UNSIGNED +#endif + } + + /** + * Base class of particles - represents a particle with position, + * an ID number and a variable number of properties. This class + * can be extended to include data related to a particle by the property + * manager. + * + * + * @ingroup Particle + * + */ + template + class Particle + { + public: + /** + * Empty constructor for Particle, creates a particle at the + * origin. + */ + Particle (); + + /** + * Constructor for Particle, creates a particle with the specified + * ID at the specified location. Note that Aspect + * does not check for duplicate particle IDs so the generator must + * make sure the IDs are unique over all processes. + * + * @param[in] new_location Initial location of particle. + * @param[in] new_reference_location Initial location of the particle + * in the coordinate system of the reference cell. + * @param[in] new_id Globally unique ID number of particle. + */ + Particle (const Point &new_location, + const Point &new_reference_location, + const types::particle_index new_id); + + /** + * Copy-Constructor for Particle, creates a particle with exactly the + * state of the input argument. Note that since each particle has a + * handle for a certain piece of the property memory, and is responsible + * for registering and freeing this memory in the property pool this + * constructor registers a new chunk, and copies the properties. + */ + Particle (const Particle &particle); + + /** + * Constructor for Particle, creates a particle from a data vector. + * This constructor is usually called after sending a particle to a + * different process. + * + * @param[in,out] begin_data A pointer to a memory location from which + * to read the information that completely describes a particle. This + * class then de-serializes its data from this memory location and + * advance the pointer accordingly. + * + * @param[in,out] new_property_pool A property pool that is used to + * allocate the property data used by this particle. + */ + Particle (const void *&begin_data, + PropertyPool &new_property_pool); + + /** + * Move constructor for Particle, creates a particle from an existing + * one by stealing its state. + */ + Particle (Particle &&particle); + + /** + * Copy assignment operator. + */ + Particle &operator=(const Particle &particle); + + /** + * Move assignment operator. + */ + Particle &operator=(Particle &&particle); + + /** + * Destructor. Releases the property handle if it is valid, and + * therefore frees that memory space for other particles. (Note: + * the memory is managed by the property_pool, and the memory is not + * deallocated by this function, it is kept in reserve for other + * particles). + */ + ~Particle (); + + /** + * Write particle data into a data array. The array is expected + * to be large enough to take the data, and the void pointer should + * point to the first element in which the data should be written. This + * function is meant for serializing all particle properties and + * afterwards de-serializing the properties by calling the appropriate + * constructor Particle(void *&data, const unsigned int data_size); + * + * @param [in,out] data The memory location to write particle data + * into. This pointer points to the begin of the memory, in which the + * data will be written and it will be advanced by the serialized size + * of this particle. + */ + void + write_data(void *&data) const; + + /** + * Set the location of this particle. Note that this does not check + * whether this is a valid location in the simulation domain. + * + * @param [in] new_loc The new location for this particle. + */ + void + set_location (const Point &new_loc); + + /** + * Get the location of this particle. + * + * @return The location of this particle. + */ + const Point & + get_location () const; + + /** + * Set the reference location of this particle. Note that this does not + * check whether this is a valid location in the simulation domain. + * + * @param [in] new_loc The new reference location for this particle. + */ + void + set_reference_location (const Point &new_loc); + + /** + * Get the reference location of this particle in its current cell. + * + * @return The reference location of this particle. + */ + const Point & + get_reference_location () const; + + /** + * Get the ID number of this particle. + * + * @return The id of this particle. + */ + types::particle_index + get_id () const; + + /** + * Tell the particle where to store its properties (even if it does not + * own properties). Usually this is only done once per particle, but + * since the particle generator does not know about the properties + * we want to do it not at construction time. Another use for this + * function is after particle transfer to a new process. + */ + void + set_property_pool(PropertyPool &property_pool); + + /** + * Set the properties of this particle. + * + * @param [in] new_properties A vector containing the + * new properties for this particle. + */ + void + set_properties (const std::vector &new_properties); + + /** + * Get write-access to properties of this particle. + * + * @return An ArrayView of the properties of this particle. + */ + const ArrayView + get_properties (); + + /** + * Get read-access to properties of this particle. + * + * @return An ArrayView of the properties of this particle. + */ + const ArrayView + get_properties () const; + + /** + * Serialize the contents of this class. + */ + template + void serialize (Archive &ar, const unsigned int version); + + private: + /** + * Current particle location. + */ + Point location; + + /** + * Current particle location in the reference cell. + * Storing this reduces the number of times we need to compute this + * location, which takes a significant amount of computing time. + */ + Point reference_location; + + /** + * Globally unique ID of particle. + */ + types::particle_index id; + + /** + * A pointer to the property pool. Necessary to translate from the + * handle to the actual memory locations. + */ + PropertyPool *property_pool; + + /** + * A handle to all particle properties + */ + PropertyPool::Handle properties; + }; + + /* -------------------------- inline and template functions ---------------------- */ + + template + template + void Particle::serialize (Archive &ar, const unsigned int) + { + ar &location + & id + & properties + ; + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif + diff --git a/include/deal.II/particle/property_pool.h b/include/deal.II/particle/property_pool.h new file mode 100644 index 0000000000..99dc52c028 --- /dev/null +++ b/include/deal.II/particle/property_pool.h @@ -0,0 +1,88 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__particle_property_pool_h +#define dealii__particle_property_pool_h + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Particle +{ + /** + * This class manages the memory space in which particles store their + * properties. Because this is dynamic memory and every particle needs the + * same amount it is more efficient to let this be handled by a central + * manager that does not need to allocate/deallocate memory every time a + * particle is constructed/destroyed. + */ + class PropertyPool + { + public: + /** + * Typedef for the handle that is returned to the particles, and that + * uniquely identifies the slot of memory that is reserved for this + * particle. + */ + typedef double *Handle; + + /** + * Define a default (invalid) value for handles. + */ + static const Handle invalid_handle; + + /** + * Constructor. Stores the number of properties per reserved slot. + */ + PropertyPool (const unsigned int n_properties_per_slot); + + /** + * Returns a new handle that allows accessing the reserved block + * of memory. + */ + Handle allocate_properties_array (); + + /** + * Mark the properties corresponding to the handle @p handle as + * deleted. After calling this function calling get_properties() with + * the freed handle causes undefined behavior. + */ + void deallocate_properties_array (const Handle handle); + + /** + * Return an ArrayView to the properties that correspond to the given + * handle @p handle. + */ + ArrayView get_properties (const Handle handle); + + /** + * Reserves the dynamic memory needed for storing the properties of + * @p size particles. + */ + void reserve(const std::size_t size); + + private: + /** + * The number of properties that are reserved per particle. + */ + const unsigned int n_properties; + }; + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 6b772be014..e21a353811 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -48,6 +48,7 @@ ADD_SUBDIRECTORY(integrators) ADD_SUBDIRECTORY(matrix_free) ADD_SUBDIRECTORY(meshworker) ADD_SUBDIRECTORY(opencascade) +ADD_SUBDIRECTORY(particle) ADD_SUBDIRECTORY(physics) ADD_SUBDIRECTORY(non_matching) ADD_SUBDIRECTORY(sundials) diff --git a/source/particle/CMakeLists.txt b/source/particle/CMakeLists.txt new file mode 100644 index 0000000000..99db585c88 --- /dev/null +++ b/source/particle/CMakeLists.txt @@ -0,0 +1,33 @@ +## --------------------------------------------------------------------- +## +## Copyright (C) 2016 by the deal.II authors +## +## This file is part of the deal.II library. +## +## The deal.II library is free software; you can use it, redistribute +## it, and/or modify it under the terms of the GNU Lesser General +## Public License as published by the Free Software Foundation; either +## version 2.1 of the License, or (at your option) any later version. +## The full text of the license can be found in the file LICENSE at +## the top level of the deal.II distribution. +## +## --------------------------------------------------------------------- + + +INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) + +SET(_src + particle.cc + property_pool.cc + ) + +SET(_inst + particle.inst.in + ) + +FILE(GLOB _header + ${CMAKE_SOURCE_DIR}/include/deal.II/particle/*.h + ) + +DEAL_II_ADD_LIBRARY(obj_particle OBJECT ${_src} ${_header} ${_inst}) +EXPAND_INSTANTIATIONS(obj_particle "${_inst}") diff --git a/source/particle/particle.cc b/source/particle/particle.cc new file mode 100644 index 0000000000..f842b41148 --- /dev/null +++ b/source/particle/particle.cc @@ -0,0 +1,268 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Particle +{ + template + Particle::Particle () + : + location (), + reference_location(), + id (0), + property_pool(NULL), + properties(PropertyPool::invalid_handle) + { + } + + + template + Particle::Particle (const Point &new_location, + const Point &new_reference_location, + const types::particle_index new_id) + : + location (new_location), + reference_location (new_reference_location), + id (new_id), + property_pool(NULL), + properties (PropertyPool::invalid_handle) + { + } + + template + Particle::Particle (const Particle &particle) + : + location (particle.get_location()), + reference_location (particle.get_reference_location()), + id (particle.get_id()), + property_pool(particle.property_pool), + properties ((property_pool != NULL) ? property_pool->allocate_properties_array() : PropertyPool::invalid_handle) + { + if (property_pool != NULL) + { + const ArrayView my_properties = property_pool->get_properties(properties); + + if (my_properties.size() != 0) + { + const ArrayView their_properties = particle.get_properties(); + + std::copy(&their_properties[0],&their_properties[0]+their_properties.size(),&my_properties[0]); + } + } + } + + + template + Particle::Particle (const void *&data, + PropertyPool &new_property_pool) + { + const types::particle_index *id_data = static_cast (data); + id = *id_data++; + const double *pdata = reinterpret_cast (id_data); + + for (unsigned int i = 0; i < dim; ++i) + location(i) = *pdata++; + + for (unsigned int i = 0; i < dim; ++i) + reference_location(i) = *pdata++; + + property_pool = &new_property_pool; + properties = property_pool->allocate_properties_array(); + + // See if there are properties to load + const ArrayView particle_properties = property_pool->get_properties(properties); + for (unsigned int i = 0; i < particle_properties.size(); ++i) + particle_properties[i] = *pdata++; + + data = static_cast (pdata); + } + +#ifdef DEAL_II_WITH_CXX11 + + template + Particle::Particle (Particle &&particle) + : + location (particle.location), + reference_location(particle.reference_location), + id (particle.id), + property_pool(particle.property_pool), + properties(particle.properties) + { + particle.properties = PropertyPool::invalid_handle; + } + + template + Particle & + Particle::operator=(const Particle &particle) + { + if (this != &particle) + { + location = particle.location; + reference_location = particle.reference_location; + id = particle.id; + property_pool = particle.property_pool; + + if (property_pool != NULL) + { + properties = property_pool->allocate_properties_array(); + const ArrayView their_properties = particle.get_properties(); + + if (their_properties.size() != 0) + { + const ArrayView my_properties = property_pool->get_properties(properties); + std::copy(&their_properties[0],&their_properties[0]+their_properties.size(),&my_properties[0]); + } + } + else + properties = PropertyPool::invalid_handle; + } + return *this; + } + + template + Particle & + Particle::operator=(Particle &&particle) + { + if (this != &particle) + { + location = particle.location; + reference_location = particle.reference_location; + id = particle.id; + property_pool = particle.property_pool; + properties = particle.properties; + particle.properties = PropertyPool::invalid_handle; + } + return *this; + } +#endif + + template + Particle::~Particle () + { + if (properties != PropertyPool::invalid_handle) + property_pool->deallocate_properties_array(properties); + } + + template + void + Particle::write_data (void *&data) const + { + types::particle_index *id_data = static_cast (data); + *id_data = id; + ++id_data; + double *pdata = reinterpret_cast (id_data); + + // Write location data + for (unsigned int i = 0; i < dim; ++i,++pdata) + *pdata = location(i); + + // Write reference location data + for (unsigned int i = 0; i < dim; ++i,++pdata) + *pdata = reference_location(i); + + // Write property data + const ArrayView particle_properties = property_pool->get_properties(properties); + for (unsigned int i = 0; i < particle_properties.size(); ++i,++pdata) + *pdata = particle_properties[i]; + + data = static_cast (pdata); + } + + template + void + Particle::set_location (const Point &new_loc) + { + location = new_loc; + } + + template + const Point & + Particle::get_location () const + { + return location; + } + + template + void + Particle::set_reference_location (const Point &new_loc) + { + reference_location = new_loc; + } + + template + const Point & + Particle::get_reference_location () const + { + return reference_location; + } + + template + types::particle_index + Particle::get_id () const + { + return id; + } + + template + void + Particle::set_property_pool (PropertyPool &new_property_pool) + { + property_pool = &new_property_pool; + } + + template + void + Particle::set_properties (const std::vector &new_properties) + { + if (properties == PropertyPool::invalid_handle) + properties = property_pool->allocate_properties_array(); + + const ArrayView old_properties = property_pool->get_properties(properties); + + std::copy(new_properties.begin(),new_properties.end(),&old_properties[0]); + } + + template + const ArrayView + Particle::get_properties () const + { + Assert(property_pool != NULL, + ExcInternalError()); + + return property_pool->get_properties(properties); + } + + + template + const ArrayView + Particle::get_properties () + { + Assert(property_pool != NULL, + ExcInternalError()); + + return property_pool->get_properties(properties); + } +} + +DEAL_II_NAMESPACE_CLOSE + +DEAL_II_NAMESPACE_OPEN + +#include "particle.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/particle/particle.inst.in b/source/particle/particle.inst.in new file mode 100644 index 0000000000..0c888d702c --- /dev/null +++ b/source/particle/particle.inst.in @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) +{ + namespace Particle + \{ + template + class Particle ; + \} +} diff --git a/source/particle/property_pool.cc b/source/particle/property_pool.cc new file mode 100644 index 0000000000..1aa1adaa05 --- /dev/null +++ b/source/particle/property_pool.cc @@ -0,0 +1,64 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Particle +{ + const PropertyPool::Handle PropertyPool::invalid_handle = NULL; + + + PropertyPool::PropertyPool (const unsigned int n_properties_per_slot) + : + n_properties (n_properties_per_slot) + {} + + + + PropertyPool::Handle + PropertyPool::allocate_properties_array () + { + return new double[n_properties]; + } + + + + void + PropertyPool::deallocate_properties_array (Handle handle) + { + delete[] handle; + } + + + + ArrayView + PropertyPool::get_properties (const Handle handle) + { + return ArrayView(handle, n_properties); + } + + + void + PropertyPool::reserve(const std::size_t size) + { + (void)size; + } +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/particle/CMakeLists.txt b/tests/particle/CMakeLists.txt new file mode 100644 index 0000000000..d33eafa20b --- /dev/null +++ b/tests/particle/CMakeLists.txt @@ -0,0 +1,4 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) +INCLUDE(../setup_testsubproject.cmake) +PROJECT(testsuite CXX) +DEAL_II_PICKUP_TESTS() diff --git a/tests/particle/particle_01.cc b/tests/particle/particle_01.cc new file mode 100644 index 0000000000..fa8a009385 --- /dev/null +++ b/tests/particle/particle_01.cc @@ -0,0 +1,51 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// check the creation and destruction of particles + +#include "../tests.h" +#include +#include +#include + + +template +void test () +{ + { + Particle::Particle particle; + + deallog << "Particle location: " << particle.get_location() << std::endl; + + Point position; + position(0) = 1.0; + particle.set_location(position); + + deallog << "Particle location: " << particle.get_location() << std::endl; + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + initlog(); + test<2>(); + test<3>(); +} diff --git a/tests/particle/particle_01.output b/tests/particle/particle_01.output new file mode 100644 index 0000000000..cbfa261dee --- /dev/null +++ b/tests/particle/particle_01.output @@ -0,0 +1,7 @@ + +DEAL::Particle location: 0.00000 0.00000 +DEAL::Particle location: 1.00000 0.00000 +DEAL::OK +DEAL::Particle location: 0.00000 0.00000 0.00000 +DEAL::Particle location: 1.00000 0.00000 0.00000 +DEAL::OK