manager.h

Go to the documentation of this file.

/*
 Copyright (C) 2012 - 2021 by the authors of the ASPECT code.

 This file is part of ASPECT.

 ASPECT is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2, or (at your option)
 any later version.

 ASPECT is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with ASPECT; see the file LICENSE.  If not see
 <http://www.gnu.org/licenses/>.
 */

#ifndef _aspect_particle_manager_h
#define _aspect_particle_manager_h

#include <aspect/global.h>

#include <deal.II/particles/particle.h>
#include <deal.II/particles/particle_accessor.h>
#include <deal.II/particles/particle_iterator.h>
#include <deal.II/particles/particle_handler.h>
#include <deal.II/particles/property_pool.h>
#include <aspect/particle/distribution.h>

#include <deal.II/matrix_free/fe_point_evaluation.h>

#include <aspect/particle/generator/interface.h>
#include <aspect/particle/integrator/interface.h>
#include <aspect/particle/interpolator/interface.h>
#include <aspect/particle/property/interface.h>

#include <aspect/simulator_access.h>
#include <aspect/simulator_signals.h>
#include <aspect/solution_evaluator.h>

#include <deal.II/base/timer.h>
#include <deal.II/base/array_view.h>

#include <boost/serialization/unique_ptr.hpp>

#include <random>
#include <vector>

namespace aspect
{
  template <int dim>
  struct SimulatorSignals;

  namespace Particle
  {
    using namespace ::Particles;

    namespace Generator
    {
      template <int dim>
      class Interface;
    }


    namespace Property
    {
      template <int dim>
      class Manager;
    }

    template <int dim>
    class Manager : public SimulatorAccess<dim>
    {
      public:
        Manager();

        ~Manager() override;

        Manager(Manager &&) noexcept;

        void initialize();

        void update();

        const Property::Manager<dim> &
        get_property_manager() const;

        const Particles::ParticleHandler<dim> &
        get_particle_handler() const;

        Particles::ParticleHandler<dim> &
        get_particle_handler();

        void copy_particle_handler (const Particles::ParticleHandler<dim> &from_particle_handler,
                                    Particles::ParticleHandler<dim> &to_particle_handler) const;

        void backup_particles ();

        void restore_particles ();


        void setup_initial_state ();

        const Interpolator::Interface<dim> &
        get_interpolator() const;

        void generate_particles();
        void initialize_particles();

        void advance_timestep();

        types::particle_index n_global_particles() const;

        void
        connect_to_signals(aspect::SimulatorSignals<dim> &signals);

#if DEAL_II_VERSION_GTE(9,6,0)
        unsigned int
        cell_weight(const typename parallel::distributed::Triangulation<dim>::cell_iterator &cell,
                    const CellStatus status);
#else
        unsigned int
        cell_weight(const typename parallel::distributed::Triangulation<dim>::cell_iterator &cell,
                    const typename parallel::distributed::Triangulation<dim>::CellStatus status);
#endif

        void update_particles();

        template <class Archive>
        void serialize (Archive &ar, const unsigned int version);

        virtual
        void
        save (std::ostringstream &os) const;

        virtual
        void
        load (std::istringstream &is);

        static
        void
        declare_parameters (ParameterHandler &prm);

        virtual
        void
        parse_parameters (ParameterHandler &prm, const unsigned int particle_manager);

      private:
        struct ParticleLoadBalancing
        {
          enum Kind
          {
            no_balancing = 0x0,
            remove_particles = 0x1,
            add_particles = 0x2,
            repartition = 0x4,
            remove_and_add_particles = remove_particles | add_particles
          };
        };

        enum class DeletionAlgorithm
        {
          random,
          point_density_function,
        };

        enum class AdditionAlgorithm
        {
          random,
          histogram,
          point_density_function,
          monte_carlo,
        };

        unsigned int addition_granularity_pdf;

        unsigned int addition_granularity_histogram;

        typename ParticlePDF<dim>::KernelFunction kernel_function;

        double bandwidth;

        std::unique_ptr<Generator::Interface<dim>> generator;

        std::unique_ptr<Integrator::Interface<dim>> integrator;

        std::mt19937 random_number_generator;

        std::unique_ptr<Interpolator::Interface<dim>> interpolator;

        std::unique_ptr<Particles::ParticleHandler<dim>> particle_handler;

        Particles::ParticleHandler<dim> particle_handler_backup;

        std::unique_ptr<Property::Manager<dim>> property_manager;

        typename ParticleLoadBalancing::Kind particle_load_balancing;

        DeletionAlgorithm deletion_algorithm;

        AdditionAlgorithm addition_algorithm;

        unsigned int min_particles_per_cell;

        unsigned int max_particles_per_cell;

        unsigned int particle_weight;

        std::map<types::subdomain_id, unsigned int>
        get_subdomain_id_to_neighbor_map() const;

        void
        apply_particle_per_cell_bounds();

        void advect_particles();

        void
        local_initialize_particles(const typename ParticleHandler<dim>::particle_iterator &begin_particle,
                                   const typename ParticleHandler<dim>::particle_iterator &end_particle);

        void
        local_update_particles(Property::ParticleUpdateInputs<dim> &inputs,
                               small_vector<Point<dim>> &positions,
                               const std::vector<EvaluationFlags::EvaluationFlags> &evaluation_flags,
                               SolutionEvaluator<dim> &evaluator);

        void
        local_advect_particles(const typename DoFHandler<dim>::active_cell_iterator &cell,
                               const typename ParticleHandler<dim>::particle_iterator &begin_particle,
                               const typename ParticleHandler<dim>::particle_iterator &end_particle,
                               SolutionEvaluator<dim> &evaluators);

        void
        connect_particle_handler_signals(aspect::SimulatorSignals<dim> &signals,
                                         ParticleHandler<dim> &particle_handler,
                                         const bool connect_to_checkpoint_signals = true) const;

        std::vector<typename Particles::ParticleHandler<dim>::particle_iterator_range>
        get_neighboring_particle_ranges(const typename Triangulation<dim>::active_cell_iterator &cell,
                                        const typename Particles::ParticleHandler<dim> &particle_handler,
                                        typename GridTools::Cache<dim> &grid_cache);
    };

    /* -------------------------- inline and template functions ---------------------- */

    template <int dim>
    template <class Archive>
    void Manager<dim>::serialize (Archive &ar, const unsigned int)
    {
      // Note that although Boost claims to handle serialization of pointers
      // correctly, at least for the case of unique_ptr it seems to not work.
      // It works correctly when archiving the content of the pointer instead.
      ar
      &(*particle_handler)
      ;
    }
  }
}

#endif