simulation.hpp

#ifndef __SIMULATIONS_UNIT_HPP__
#define __SIMULATIONS_UNIT_HPP__
 
#include <Kokkos_Core.hpp>
#include <Kokkos_ScatterView.hpp>
#include <biocma_cst_config.hpp>
#include <cassert>
#include <cma_utils/alias.hpp>
#include <common/common.hpp>
#include <common/logger.hpp>
#include <cstddef>
#include <impl/Kokkos_Half_FloatingPointWrapper.hpp>
#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;
 
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;
    void update_feed(double t, double d_t, bool update_scalar = true) noexcept;
    // void update(CmaUtils::IterationState&& new_state);
    void updateHydro(const CmaUtils::IterationStatePtrType& newstate);
 
    void updateMCHydro(std::span<const double> vl,
                       std::span<const std::size_t> neighors_flat,
                       std::span<const double> proba_flat,
                       std::span<const double> out_flows);
 
    bool checkScalar() const;
 
    // Memory management
    void reset();
 
    template <typename Space, ModelType Model>
    KernelInline::CycleFunctors<Space, Model>
    init_functors(MC::ParticlesContainer<Model> container);
 
  private:
    void updateScalarHydro(const CmaUtils::IterationStatePtrType& newstate);
 
    void setVolumes(const CmaUtils::IterationStatePtrType& newstate) const;
    void setLiquidFlow(CmaUtils::PreCalculatedHydroState* _flows_l);
    void setGasFlow(CmaUtils::PreCalculatedHydroState* _flows_g);
    void post_init_compartments();
 
    [[nodiscard]] MC::KernelConcentrationType getkernel_concentration() const;
 
    MC::ContributionView contribs_scatter;
    MapProbes probes;
 
    // Attributes
    // 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
    bool f_reaction = true; // FIXME
    void scatter_contribute();
    void set_kernel_contribs_to_host();
 
    [[nodiscard]] kernelContribution get_kernel_contribution() const;
 
    void post_init_concentration(const ScalarInitializer& scalar_init);
 
    // std::unique_ptr<ScalarSimulationProxyBase> scalar_proxy;
 
    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 <ModelType Model>
    void post_cycle(MC::ParticlesContainer<Model>& container,
                    auto& cycle_functors);
 
    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& cycle_functors)
  {
    PROFILE_SECTION("Simulation::pre_cycle")
    // 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();
 
    cycle_functors.update(d_t, container, this->move_info);
  }
 
  template <typename Space, ModelType Model>
  KernelInline::CycleFunctors<Space, Model>
  SimulationUnit::init_functors(MC::ParticlesContainer<Model> container)
  {
 
    return KernelInline::CycleFunctors<Space, Model>(
        container,
        mc_unit->rng.random_pool,
        getkernel_concentration(),
        contribs_scatter,
        mc_unit->events,
        move_info,
        probes[ProbeType ::LeavingTime]);
  }
 
  void SimulationUnit::cycleProcess(auto& container,
                                    double d_t,
                                    auto& cycle_functors)
  {
    PROFILE_SECTION("cycleProcess")
    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, cycle_functors);
 
    if (f_reaction)
    {
      this->contribs_scatter.reset();
      cycle_functors.launch_model(n_particle);
    }
 
    cycle_functors.move_space.fence();
    if (cycle_functors.move_kernel.need_launch())
    {
 
      cycle_functors.launch_move(n_particle);
    }
 
    post_cycle<CurrentModel>(container, cycle_functors);
  }
 
  template <ModelType Model>
  void SimulationUnit::post_cycle(MC::ParticlesContainer<Model>& container,
                                  auto& cycle_functors)
  {
    PROFILE_SECTION("Simulation::post_cycle")
    Kokkos::fence();
    this->scatter_contribute();
    auto [host_red, host_out_counter] = cycle_functors.get_host_reduction();
 
    container.update_and_remove_inactive(host_out_counter, host_red.dead_total);
 
    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__