simulation.hpp

#ifndef __SIMULATIONS_UNIT_HPP__
#define __SIMULATIONS_UNIT_HPP__
 
#include "simulation/kernels/model_kernel.hpp"
#include <Kokkos_Core.hpp>
#include <Kokkos_ScatterView.hpp>
#include <biocma_cst_config.hpp>
#include <cassert>
#include <cma_utils/iteration_state.hpp>
#include <common/common.hpp>
#include <common/logger.hpp>
#include <cstddef>
#include <cstdint>
#include <mc/domain.hpp>
#include <mc/events.hpp>
#include <mc/prng/prng.hpp>
#include <mc/unit.hpp>
#include <memory>
#include <optional>
#include <simulation/alias.hpp>
#include <simulation/feed_descriptor.hpp>
#include <simulation/kernels/kernels.hpp>
#include <simulation/mass_transfer.hpp>
#include <simulation/probe.hpp>
#include <simulation/scalar_initializer.hpp>
 
// TODO Clean
static constexpr size_t trigger_const_particle_number = 1e6;
static constexpr bool f_reaction = true; // FIXME
 
namespace CmaUtils
{
  class PreCalculatedHydroState;
} // namespace CmaUtils

namespace Simulation
{
 
  struct Dimensions
  {
    std::size_t n_species{};
    std::size_t n_compartment{};
 
    template <class Archive> void serialize(Archive& archive)
    {
      archive(n_species, n_compartment);
    }
  };
 
  class ScalarSimulation;
 
  class SimulationUnit
  {
  public:
    SimulationUnit(std::unique_ptr<MC::MonteCarloUnit>&& _unit,
                   const ScalarInitializer& scalar_init,
                   std::optional<Feed::SimulationFeed> _feed = std::nullopt);
 
    ~SimulationUnit();
 
    SimulationUnit(SimulationUnit&& other) noexcept;
    SimulationUnit(const SimulationUnit& other) = delete;
    SimulationUnit& operator=(SimulationUnit&& rhs) = delete;
    SimulationUnit& operator=(const SimulationUnit& rhs) = delete;
 
    std::unique_ptr<MC::MonteCarloUnit> mc_unit; // NOLINT
 
    [[nodiscard]] const Simulation::Feed::SimulationFeed& get_feed() const;
    ProbeAutogeneratedBuffer& get_probes();
 
    // std::span<double> getContributionData_mut();
 
    // Setters
    void setProbes(ProbeAutogeneratedBuffer&& _probes);
    void setMtrModel(MassTransfer::Type::MtrTypeVariant&& variant);
    void setLogger(std::shared_ptr<IO::Logger>);
 
    // Getters
    [[nodiscard]] double& get_start_time_mut();
    [[nodiscard]] double& get_end_time_mut();
    [[nodiscard]] const CmaUtils::IterationState& getState() const;
    [[nodiscard]] Dimensions getDimensions() const noexcept;
    [[nodiscard]] std::span<double> getCliqData() const;
    [[nodiscard]] std::optional<std::span<const double>> getCgasData() const;
    [[nodiscard]] std::span<const double> getContributionData() const;
    [[nodiscard]] bool two_phase_flow() const;
    [[nodiscard]] std::optional<std::span<const double>> getMTRData() const;
    [[nodiscard]] std::size_t dead_counter() const;
    std::span<double> getContributionData_mut();
 
    // Simulation methods
 
    void cycleProcess(auto&& container, double d_t, auto& _functors);
    void step(double d_t) const;
    void reduceContribs(std::span<const double> data, size_t n_rank) const;
 
    //[[deprecated("perf:not useful")]] void
    // reduceContribs_per_rank(std::span<const double> data) const;
    //
    void clearContribution() const noexcept;
    auto init_functor(double d_t, auto&& container);
    void update_feed(double t, double d_t, bool update_scalar = true) noexcept;
    void update(CmaUtils::IterationState&& new_state);
 
    // Memory management
    void clear_mc();
    void reset();
    void post_init_compartments();
 
    template <typename Space, ModelType Model>
    KernelInline::Functors<Space, Model>
    init_functors(MC::ParticlesContainer<Model> container);
 
  private:
    void setVolumes() const;
    void setLiquidFlow(CmaUtils::PreCalculatedHydroState* _flows_l);
    void setGasFlow(CmaUtils::PreCalculatedHydroState* _flows_g);
 
    [[nodiscard]] MC::KernelConcentrationType getkernel_concentration() const;
 
    MC::ContributionView contribs_scatter;
    MapProbes probes;
 
    // Attributes
    std::size_t internal_counter_dead;
    uint64_t minimum_dead_particle_removal;
    // ProbeAutogeneratedBuffer probes;
    Simulation::Feed::SimulationFeed feed;
 
    CmaUtils::IterationState state;
    bool const_number_simulation = true;
    bool is_two_phase_flow;
    double starting_time =
        0.;            // Not used within calculation, only for export purposes
    double end_time{}; // Not used within calculation, only for export purposes
 
    void scatter_contribute();
    void set_kernel_contribs_to_host();
    // Bounce methods to pimpl
    [[nodiscard]] DiagonalView<ComputeSpace> get_kernel_diagonal() const;
    [[nodiscard]] CumulativeProbabilityView<ComputeSpace>
    get_kernel_cumulative_proba() const;
    [[nodiscard]] kernelContribution get_kernel_contribution() const;
 
    void post_init_concentration(const ScalarInitializer& scalar_init);
    void post_init_concentration_functor(const ScalarInitializer& scalar_init);
    void post_init_concentration_file(const ScalarInitializer& scalar_init);
 
    std::shared_ptr<ScalarSimulation> liquid_scalar;
    std::shared_ptr<ScalarSimulation> gas_scalar;
    MassTransfer::MassTransferModel
        mt_model; // TODO add default null value (no model)
    KernelInline::MoveInfo<ComputeSpace> move_info;
 
    std::shared_ptr<IO::Logger> logger;
 
    template <typename Space, ModelType Model>
    void post_cycle(
        MC::ParticlesContainer<Model>& container,
        std::size_t n_particle,
        KernelInline::Functors<Space, Model>::move_reducer_view_type out_total,
        KernelInline::Functors<Space, Model>::cycle_reducer_view_type
            reducer_type);
 
    template <ModelType Model>
    void pre_cycle(MC::ParticlesContainer<Model>& container,
                   double d_t,
                   auto& functors);
  };
 
  template <ModelType Model>
  void SimulationUnit::pre_cycle(MC::ParticlesContainer<Model>& container,
                                 double d_t,
                                 auto& functors)
  {
    this->contribs_scatter.reset();
    this->move_info.cumulative_probability = get_kernel_cumulative_proba();
    this->move_info.diag_transition = get_kernel_diagonal();
    this->move_info.neighbors = mc_unit->domain.getNeighbors();
 
    //TODO map through flowmap and find max flow, then true if max>0
    bool enable_move = move_info.liquid_volume.size() > 1;
    bool enable_leave = move_info.leaving_flow.size() != 0;
 
    // FIXME: cycle_kernel: need to update container because we change:
    // n_used_element counter which is size_t outside of kernel. As functor owns
    // a copy of container this counter is not update between iterations.
    // 1. Use n_used_element as a reference counter type
    // 2. Manually update counter in update function (dirty way)
 
    functors.cycle_kernel.update(d_t, container);
 
    // TODO: Why need to update all views (where did we lost the refcount ? )
    functors.move_kernel.update(d_t,
                                container.n_particles(),
                                this->move_info,
                                container.position,
                                container.status,
                                container.ages,
                                enable_move,
                                enable_leave);
  }
 
  template <typename Space, ModelType Model>
  KernelInline::Functors<Space, Model>
  SimulationUnit::init_functors(MC::ParticlesContainer<Model> container)
  {
    auto local_rng = mc_unit->rng;
    auto events = mc_unit->events;
    auto reaction_functor =
        Simulation::KernelInline::CycleFunctor<Model>(container,
                                                      local_rng.random_pool,
                                                      getkernel_concentration(),
                                                      contribs_scatter,
                                                      events);
    auto move_functor =
        Simulation ::KernelInline ::MoveFunctor(container.position,
                                                container.status,
                                                move_info,
                                                local_rng.random_pool,
                                                events,
                                                probes[ProbeType ::LeavingTime],
                                                container.ages);
 
    return KernelInline::Functors<Space, Model>(std::move(reaction_functor),
                                                std::move(move_functor));
  }
 
  void
  SimulationUnit::cycleProcess(auto&& container, double d_t, auto& functors)
  {
 
    PROFILE_SECTION("cycleProcess")
 
    using cycle_functor_type = decltype(functors.cycle_kernel);
    using CurrentModel =
        typename std::remove_reference<decltype(container)>::type::UsedModel;
    const size_t n_particle = container.n_particles();
    if (n_particle == 0)
    {
      return;
    }
 
    pre_cycle(container, d_t, functors);
 
    if (functors.move_kernel.need_launch())
    {
      constexpr bool is_reduce = true;
      const auto _policy =
          MC::get_policy(functors.move_kernel, n_particle, is_reduce);
      Kokkos ::parallel_reduce(
          "cycle_move", _policy, functors.move_kernel, functors.move_reducer);
    }
 
    if (f_reaction)
    {
      constexpr bool is_reduce = true;
      const auto _policy =
          MC::get_policy<cycle_functor_type, KernelInline::TagSecondPass>(
              functors.cycle_kernel, n_particle, is_reduce);
      const auto scatter_policy =
          MC::get_policy<cycle_functor_type, KernelInline::TagFirstPass>(
              functors.cycle_kernel, n_particle, false);
 
      Kokkos::parallel_reduce(
          "cycle_model",
          _policy,
          functors.cycle_kernel,
          KernelInline::CycleReducer<ComputeSpace>(functors.cycle_reducer));
      Kokkos::fence(); // TODO needed ?
      Kokkos::parallel_for(
          "cycle_scatter", scatter_policy, functors.cycle_kernel);
    }
 
    post_cycle<ComputeSpace, CurrentModel>(
        container, n_particle, functors.move_reducer, functors.cycle_reducer);
  }
 
  template <typename Space, ModelType Model>
  void SimulationUnit::post_cycle(
      MC::ParticlesContainer<Model>& container,
      std::size_t n_particle,
      KernelInline::Functors<Space, Model>::move_reducer_view_type out_total,
      KernelInline::Functors<Space, Model>::cycle_reducer_view_type
          reducer_type)
  {
    const auto threshold = std::max(
        minimum_dead_particle_removal,
        static_cast<uint64_t>(static_cast<double>(n_particle) *
                              AutoGenerated::dead_particle_ratio_threshold));
 
    Kokkos::fence();
 
    this->scatter_contribute();
 
    const auto host_red =
        Kokkos::create_mirror_view_and_copy(HostSpace(), reducer_type)();
 
    const auto host_out_counter =
        Kokkos::create_mirror_view_and_copy(HostSpace(), out_total)();
 
    internal_counter_dead = container.update_and_clean_dead(
        host_out_counter, host_red.dead_total, threshold);
 
    container.merge_buffer();
 
    if (host_red.waiting_allocation_particle != 0)
    {
      if (logger)
      {
        logger->alert("Simulation",
                      "Overflow of particle not implemented yet (ignore "
                      "_waiting_allocation_particle)");
      }
    }
 
    set_kernel_contribs_to_host();
  }
 
} // namespace Simulation
 
#endif //__SIMULATIONS_UNIT_HPP__