Simulation.h

#ifndef __SIMULATION_SIMULATION_H__
#define __SIMULATION_SIMULATION_H__

#include "../CellList/CellListFast.h"
#include "../EventCalendar/EventCalendar.h"
#include "../EventCalendar/EventType.h"
#include "../Particle/Particle.h"
#include "../RPK/RPK.h"
#include "../Vector/Vector.hpp"

#include "./CollisionRecord.h"
#include "./Statistics.h"

#include <math.h>
#include <vector>

namespace MDSimulation
{
    class AllParticles
    {
    public:

        explicit AllParticles(const int)
        { }

        inline bool operator()(const int)
        {
            return true;
        }
    };

    class LessThan
    {
    private:
        const int Upper;    // Strict upper bound of particles to accept.

    public:
        explicit LessThan(const int i) : Upper(i)
        { }

        inline bool operator()(const int pa)
        {
            return pa < Upper;
        }
    };

    class Skipping
    {
    private:
        const int Skip;     // Index value to skip.

    public:
        explicit Skipping(const int s) : Skip(s)
        { }

        inline bool operator()(const int pa)
        {
            return pa != Skip;
        }
    };

    class Simulation
    {
    private:

        // The number of radial bins to compute statistic for when writing a snapshot.
        static const int MAX_BINS = 48;

        // The number of collisions to be recorded.
        static const size_t COLLISION_RECORDS = 100000;

        typedef std::unique_ptr<EventCalendar> CalendarPtr;

        InfoStruct Info;            
        CellList Partitions;  
        CalendarPtr Calendar; 
        std::vector<Particle> Particles;   
        std::vector<double>   ConeTime;    
        std::vector<double>   WallTime;    
        double                CurrentTime; 
        int ParticleCount; 
        int MaxParticles;  
        int SnapshotCount;                 
        DoubleVector Extent;    
        DoubleVector Offset;    
        IntVector    Cells;     
        IntVector    InitUCell; 
        RPK        RPKEquation;     
        double     LastRebinRadius; 
        double     LastVelocity;    
        double     LastSnapshotTime; 
        bool       Rebounded;       
        bool       Save;        
        Statistics Stats;

        std::string OutDir;

        clock_t PerformanceClock;

        template <class Controller>
        inline void PredictEvents(const int na, const int nb, const EventType ev);

        int CellCrossingTimes(const int na, DoubleVector& times);

        double ConeCollisionTime(const int na);

        double WallCollisionTime(const int na);

        template <class Controller>
        void ScheduleParticleCollisions(const int na, Controller& con,
                                        const IntVector& low,
                                        const IntVector& high);

        template <class Controller>
        void ScheduleParticleCollisions(const int na, Controller& con);

        double ComputeIntersectionTime(const Particle& pa, const Particle& pb);

        void ProcessCellCrossing();

        void ProcessBubbleWall();

        void ProcessCollision();

        void DissociateParticle(const int na);

        void ProcessCone();

        void ProcessUpdateSystem();

        void ProcessReadPGauge();

        double Diameter2(const double vv, const Particle& p1, const Particle& p2);

        void BinSystem();

        void StartRun();

        int ComputeParticlesInSystem() const;

        void InitializeParticles();


        void PopulateCalendar();

        void WriteSnapshot();


        void WriteCollisionEnergy();

    public:

        friend std::ostream& operator<<(std::ostream& stream, const Simulation& sim);

        explicit Simulation(const InfoStruct& info, const std::string& outDir = "", const bool save = true);

        void PurgeEvents();

        bool HasRebounded() const;

        EventType NextStep();

        std::vector<Particle>::const_iterator begin();

        std::vector<Particle>::const_iterator end();

        double GetCurrentTime() const;

        unsigned long long GetCollisionCount() const;

        unsigned long long GetCellCrossingCount() const;

        unsigned long long GetWallCollisionCount() const;

        unsigned long long GetConeBoundaryCount() const;

        template <typename Writer>
        void Serialize(Writer& writer) const;

        int GetParticleCount() const;

        // void Read(const boost::property_tree::ptree& tree);
    };

    template <typename Controller>
    void Simulation::PredictEvents(const int na, const int nb, const EventType ev)
    {
        Controller con(nb);
        Particle& pa = Particles[na];
        pa.UpdateParticle(this->CurrentTime);

        DoubleVector cellTime(0.0, 0.0, 0.0);

        const int evCode  = this->CellCrossingTimes(na, cellTime);
        cellTime[evCode] += this->CurrentTime;

        // Do not repredict cone wall and bubble wall events after a
        // cell crossing event. The previous values are not invalidated by
        // these events.
        if (ev != CellCrossingEvent)
        {
            this->ConeTime[na] = this->ConeCollisionTime(na) + this->CurrentTime;

            // Do not repredict Bubble wall crossings after cone
            // events, as cone wall events, supposedly, do not alter the
            // radial velocity of the particles
            if (ev != ConeWallEvent)
            {
                this->WallTime[na] = this->WallCollisionTime(na);
            }
        }

        // Find the nearest event for this particle and schedule.
        if (cellTime[evCode] < this->ConeTime[na])
        {
            if (cellTime[evCode] < this->WallTime[na])
            {
                // Schedule a cell crossing event.
                this->Calendar->ScheduleEvent(
                        cellTime[evCode], CellCrossingEvent, na, evCode);
            }
            else
            {
                // Schedule a bubble wall event. The second particle in bubble
                // wall events is ignored.
                this->Calendar->ScheduleEvent(
                        this->WallTime[na], BubbleWallEvent, na, -1);
            }
        }
        else
        {
            if (this->ConeTime[na] < this->WallTime[na])
            {
                // Schedule a cone wall event. The second particle in cone
                // wall events is ignored.
                this->Calendar->ScheduleEvent(
                        this->ConeTime[na], ConeWallEvent, na, -1);
            }
            else
            {
                // Schedule a bubble wall event. The second particle in bubble
                // wall events is ignored.
                this->Calendar->ScheduleEvent(
                        this->WallTime[na], BubbleWallEvent, na, -1);
            }
        }


        // If the last event was a cell crossing, we can limit the search to
        // those in the direction we are heading, as we will have retained
        // the collisions from particles within the current and previous
        // cell.
        IntVector low(-1, -1, -1);
        IntVector high(1, 1, 1);

        // If a cell crossing just occurred, we can optimize by providing
        // more restricted bounds on the neighborhood search.
        if (ev == CellCrossingEvent)
        {
            if (pa.Velocity[nb] >= 0.0)
            {
                low[nb] = 1;
            }
            else
            {
                high[nb] = -1;
            }
        }

        this->ScheduleParticleCollisions<Controller>(na, con, low, high);
    }

    template <typename Controller>
    void Simulation::ScheduleParticleCollisions(const int na, Controller& con)
    {
        static const IntVector low(-1, -1, -1);
        static const IntVector high(1, 1, 1);
        this->ScheduleParticleCollisions<Controller>(na, con, low, high);
    }

    /*
     * Alright, this seems like it deserves some explanation, as it is somewhat of
     * a complex design.
     *
     * This is to cheat the need for having multiple versions of this function or
     * having a convoluted data encoding in the parameters/global variables.
     * Controller is a function object with a function call operator
     * that acts as a predicate function.
     *
     * If the predicate tests true, then the given particle is tested for
     * collision with the particle designated by <code>na</code>.
     *
     * In this sense, the parameter <code>con</code> acts as a closure containing
     * all the data needed to determine whether or not to use this particle.
     */
    template <typename Controller>
    void Simulation::ScheduleParticleCollisions(
            const int na, Controller& con,
            const IntVector& low,
            const IntVector& high)
    {
        assert(na >= 0);
        assert((unsigned) na < Particles.size());

        Particle& pa = Particles[na];

        CellList::NeighborCellIterator it = Partitions.BeginNeighborCell(na, low, high);

        for (; it.HasNext(); ++it)
        {
            if (*it != na && con(*it))
            {
                Particle& pb = Particles[*it];
                pb.UpdateParticle(this->CurrentTime);

                const double time = ComputeIntersectionTime(pa, pb);

                if (time != Constants::NEVER)
                {
                    this->Calendar->ScheduleEvent(
                            this->CurrentTime + time, CollisionEvent, na, *it);
                }
            }
        }
    }

    template <typename Writer>
    void Simulation::Serialize(Writer& writer) const
    {
        writer.String("Info")        , Info.Serialize(writer);
        writer.String("RPKEquation") , RPKEquation.Serialize(writer);

        writer.String("CurrentTime")         , writer.Double(CurrentTime);
        writer.String("TotalCollisions")     , writer.UInt64(Stats.TotalCollisions);
        writer.String("TotalCellCrossings")  , writer.UInt64(Stats.TotalCellCrossings);
        writer.String("TotalWallCollisions") , writer.UInt64(Stats.TotalWallCollisions);

        writer.String("TotalConeBoundaryCollisions");
        writer.UInt64(Stats.TotalConeBoundaryCollisions);

        writer.String("FusionRate"), writer.UInt64(Stats.FusionRate);

        // Serialize the array of particles to the document
        writer.String("Particles");
        writer.StartArray();
        std::for_each(Particles.begin(), Particles.end(), [&writer](const Particle& p)
                {
                    p.Serialize(writer);
                });
        writer.EndArray();
    }
}

#endif /* __SIMULATION_H__ */