world.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_world_h
#define _aspect_particle_world_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 <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 <deal.II/base/timer.h>
#include <deal.II/base/array_view.h>

#include <boost/serialization/unique_ptr.hpp>

namespace aspect
{
  template <int dim>
  struct SimulatorSignals;

  namespace Particle
  {
    using namespace dealii;
    using namespace ::Particles;

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


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

    namespace internal
    {
      template <int dim>
      class SolutionEvaluators
      {
        public:
          virtual ~SolutionEvaluators() = default;

          virtual void
          reinit(const typename DoFHandler<dim>::active_cell_iterator &cell,
                 const ArrayView<Point<dim>> &positions,
                 const ArrayView<double> &solution_values,
                 const UpdateFlags update_flags) = 0;

          virtual void get_solution(const unsigned int evaluation_point,
                                    Vector<double> &solution) = 0;

          virtual void get_gradients(const unsigned int evaluation_point,
                                     std::vector<Tensor<1, dim>> &gradients) = 0;

          virtual FEPointEvaluation<dim, dim> &
          get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity) = 0;

          virtual NonMatching::MappingInfo<dim> &
          get_mapping_info() = 0;
      };

      template <int dim>
      std::unique_ptr<internal::SolutionEvaluators<dim>>
      construct_solution_evaluators(const SimulatorAccess<dim> &simulator_access,
                                    const UpdateFlags update_flags);
    }

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

        ~World() override;

        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 world_index);

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

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

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

        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;

        unsigned int min_particles_per_cell;

        unsigned int max_particles_per_cell;

        unsigned int particle_weight;

        bool update_ghost_particles;

        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(const typename DoFHandler<dim>::active_cell_iterator &cell,
                               internal::SolutionEvaluators<dim> &evaluators);

        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,
                               internal::SolutionEvaluators<dim> &evaluators);

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

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

    template <int dim>
    template <class Archive>
    void World<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