/* PBHVlasov example (pbhgr project).  Weak profile-referenced 1+log lapse (devlog 29–30 Sep):
 *   d_t alpha = advec - eps alpha^p (K - fref(x) K_far(t) - 2 Theta),
 * fref = K_0(x)/K_0(far) carried as a static state variable (c_fref), K_far = far-field (FLRW) K set by the
 * level every coarse step (static member).  Shift: integrated Gamma-driver as in IntegratedMovingPunctureGauge. */
#ifndef PBHGAUGETECLYN_HPP_
#define PBHGAUGETECLYN_HPP_

#include "CCZ4Vars.hpp"
#include "DimensionDefinitions.hpp"
#include "IntegratedMovingPunctureGauge.hpp"
#include "StateVariables.hpp"

template <class deriv_t = FourthOrderDerivatives> class PBHGaugeTeclyn : public IntegratedMovingPunctureGauge<deriv_t>
{
  public:
    using base_t = IntegratedMovingPunctureGauge<deriv_t>;
    using params_t = typename base_t::params_t;
    PBHGaugeTeclyn(amrex::Real a_dx, amrex::Real a_K_far) : base_t(a_dx), m_K_far(a_K_far) {}
    //! the Gamma-driver's gauge speed is sqrt(shift_Gamma_coeff) in coordinate units whatever chi is; with the
    //! expansion-aware time step (coordinate light speed sqrt(chi_far)) the coefficient is scaled by chi_far
    void scale_shift_Gamma(amrex::Real a_factor) { this->m_params.shift_Gamma_coeff *= a_factor; }

    AMREX_GPU_DEVICE AMREX_FORCE_INLINE void calculate_rhs(int ix, int iy, int iz, const amrex::Array4<amrex::Real> &rhs,
                                                           const amrex::Array4<const amrex::Real> &state) const
    {
        const amrex::CellData<amrex::Real> &rhs_cell_data = rhs.cellData(ix, iy, iz);
        const amrex::CellData<const amrex::Real> &state_cell_data = state.cellData(ix, iy, iz);
        const CCZ4Vars vars(state_cell_data);
        const Tensor::Rank1 shift_vector({vars.shift(0), vars.shift(1), vars.shift(2)});
        const amrex::Real advec_lapse = this->m_deriv.advec_scalar(ix, iy, iz, state, shift_vector, c_lapse);
        const Tensor::Rank1 advec_shift = this->m_deriv.advec_vector(ix, iy, iz, state, shift_vector, c_shift1);
        amrex::Real eta_of_x{};
        this->compute_eta(eta_of_x, ix, iy, iz);
        const amrex::Real fref = state_cell_data[c_fref];
        rhs_cell_data[c_lapse] = this->m_params.lapse_advec_coeff * advec_lapse -
                                 this->m_params.lapse_coeff * pow(vars.lapse(), this->m_params.lapse_power) *
                                     (vars.K() - fref * m_K_far - 2.0 * vars.Theta());
        FOR (i)
        {
            rhs_cell_data[c_shift1 + i] = this->m_params.shift_advec_coeff * advec_shift(i) +
                                          this->m_params.shift_Gamma_coeff * vars.Gamma(i) - eta_of_x * vars.shift(i) -
                                          vars.B(i);
            rhs_cell_data[c_B1 + i] = 0.0;
        }
        rhs_cell_data[c_fref] = 0.0;
    }

  private:
    amrex::Real m_K_far;
};

#endif /* PBHGAUGETECLYN_HPP_ */