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 <mc/domain.hpp>
#include <mc/events.hpp>
#include <mc/prng/prng.hpp>
#include <mc/unit.hpp>
#include <memory>
#include <optional>
#include <ranges>
#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(ProbeType type);
 
    // std::span<double> getContributionData_mut();
 
    // Setters
    void setProbes(ProbeType type, 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]] 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;
    std::span<double> getContributionData_mut();
 
    // Simulation methods
    void cycleProcess(auto& container, double d_t, auto& _functors);
    void step(double d_t) const;
    void clearContribution() const noexcept;
    void update_feed(double t, double d_t, bool update_scalar = true) noexcept;
    void updateHydro(const CmaUtils::IterationStatePtrType& newstate);
    void updateMCHydro(std::span<const double> newliquid_volume,
                       std::span<const std::size_t> neighors_flat,
                       std::span<const double> proba_flat,
                       std::span<const double> out_flows) const;
 
    bool checkScalar() const;
 
    // Memory management
    void reset();
 
    template <typename Space, ModelType Model>
    KernelInline::CycleFunctors<Space, Model>
    init_functors(MC::ParticlesContainer<Model> container);
 
    std::ranges::subrange<ProbeIterator, ProbeIterator> it_probes() const;
 
  private:
    void updateScalarHydro(const CmaUtils::IterationStatePtrType& newstate);
 
    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
    Simulation::Feed::SimulationFeed feed;
 
    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)
 
    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")
 
    cycle_functors.update(
        d_t, container, this->mc_unit->domain.get_const_inner());
  }
 
  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,
        mc_unit->domain.get_const_inner(),
        probes[ProbeType ::LeavingTime],
        probes[ProbeType ::DivisionTime]);
  }
 
  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__