/* PBHVlasov example for GRTeclyn (pbhgr project, stage 2): BSSN/CCZ4 + Vlasov particles. */
#ifndef PBHVLASOVLEVEL_HPP_
#define PBHVLASOVLEVEL_HPP_

#include <memory>
#include "DefaultLevelBld.hpp"
#include "FourthOrderDerivatives.hpp"
#include "GRAmrLevel.hpp"
#include "PBHVlasovAmr.hpp"
#include "ParticleMatter.hpp"

class PBHVlasovLevel : public GRAmrLevel
{
  public:
    using GRAmrLevel::GRAmrLevel;
    template <class deriv_t = FourthOrderDerivatives> using Matter = ParticleMatter<deriv_t>;

    PBHVlasovAmr *get_vlasov_amr_ptr();
    static void variableSetUp();

    void initData() override;
    void specific_post_init() override;
    void specific_advance() override;
    //! pbh_vlasov.strang = 1: half particle step and deposit before the grid update (sources at mid-step), the second
    //! half in specific_advance(); second order in time when the sources change by per cents per step (early start)
    amrex::Real advance(amrex::Real time, amrex::Real dt, int iteration, int ncycle) override;
    void specific_eval_rhs(amrex::MultiFab &a_soln, amrex::MultiFab &a_rhs, amrex::Real a_time) override;
    void specific_update_ode(amrex::MultiFab &a_soln) override;
    void specific_post_timestep() override;
    void specific_post_regrid(int a_lbase, int a_new_finest) override;
    void specific_post_checkpoint(const std::string &a_dir, std::ostream &os) override;
    void specific_post_restart() override;
    void tag_cells(amrex::TagBoxArray &a_tag_box_array, amrex::Real a_regrid_threshold) final;
    //! expansion-aware time step with adaptive subcycling (pbh_vlasov.dt_mode = 1); otherwise GRTeclyn's fixed step
    void computeInitialDt(int finest_level, int sub_cycle, amrex::Vector<int> &n_cycle,
                          const amrex::Vector<amrex::IntVect> &ref_ratio, amrex::Vector<amrex::Real> &dt_level,
                          amrex::Real stop_time) override;
    void computeNewDt(int finest_level, int sub_cycle, amrex::Vector<int> &n_cycle,
                      const amrex::Vector<amrex::IntVect> &ref_ratio, amrex::Vector<amrex::Real> &dt_min,
                      amrex::Vector<amrex::Real> &dt_level, amrex::Real stop_time, int post_regrid_flag) override;
    //! `yoo` mode: energy and momentum density from the vacuum constraints of this level's data, then the level's
    //! particles from them (called for every level once the lattice exists)
    void init_particles_from_constraints();
    //! `yoo` mode: correction field of this level from the deposited coordinate sources and the constraint targets
    //! (mass ratio D_target / D_dep and du_i = (S_target - S_dep)_i / D_target; ghost cells from the coarser level and
    //! the neighbours); returns max |ratio - 1|.  apply_deposit_correction() applies it to the level's particles.
    amrex::Real build_deposit_correction();
    void apply_deposit_correction();
    void release_deposit_targets() { m_yoo_target.reset(); m_yoo_corr.reset(); }

  private:
    void deposit_particles();
    void set_time_steps(int finest_level, amrex::Vector<int> &n_cycle, amrex::Vector<amrex::Real> &dt_level);
    void load_table();
    // tabulated spherical initial data (isotropic radius): chi, K, A^r_r, rho, u_r
    std::vector<amrex::Real> m_tab_r, m_tab_chi, m_tab_K, m_tab_Arr, m_tab_rho, m_tab_ur;
    std::vector<amrex::Real> m_tab_sr, m_tab_st; // optional: sigma_r, sigma_t ratios (warm start)
    // yoo mode, initialisation only: E, J_i, chi, h~_ij from the constraints (1 ghost) and the correction field
    std::unique_ptr<amrex::MultiFab> m_yoo_target, m_yoo_corr;
    bool m_table_loaded{false};
};

#endif /* PBHVLASOVLEVEL_HPP_ */