/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2019 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
#ifndef SWIFT_GIZMO_MFM_HYDRO_FLUX_H
#define SWIFT_GIZMO_MFM_HYDRO_FLUX_H

#include "riemann.h"

/**
 * @brief Reset the hydrodynamical fluxes for the given particle.
 *
 * @param p Particle.
 */
__attribute__((always_inline)) INLINE static void hydro_part_reset_hydro_fluxes(
    struct part* restrict p) {

  p->flux.momentum[0] = 0.0f;
  p->flux.momentum[1] = 0.0f;
  p->flux.momentum[2] = 0.0f;
  p->flux.energy = 0.0f;
}

/**
 * @brief Reset the gravity fluxes for the given particle.
 *
 * Does nothing, since MFM does not have gravity fluxes.
 *
 * @param p Particle.
 */
__attribute__((always_inline)) INLINE static void
hydro_part_reset_gravity_fluxes(struct part* restrict p) {}

/**
 * @brief Get the fluxes for the given particle.
 *
 * @param p Particle.
 * @param flux Fluxes for the particle (array of size 5 or more).
 */
__attribute__((always_inline)) INLINE static void hydro_part_get_fluxes(
    const struct part* restrict p, float* flux) {

  flux[1] = p->flux.momentum[0];
  flux[2] = p->flux.momentum[1];
  flux[3] = p->flux.momentum[2];
  flux[4] = p->flux.energy;
}

/**
 * @brief Compute the flux for the Riemann problem with the given left and right
 * state, and interface normal, surface area and velocity.
 *
 * @param WL Left state variables.
 * @param WR Right state variables.
 * @param n_unit Unit vector of the interface.
 * @param vLR Velocity of the interface.
 * @param Anorm Surface area of the interface.
 * @param fluxes Array to store the result in (of size 5 or more).
 */
__attribute__((always_inline)) INLINE static void hydro_compute_flux(
    const float* WL, const float* WR, const float* n_unit, const float* vLR,
    const float Anorm, float* fluxes) {

  riemann_solve_for_middle_state_flux(WL, WR, n_unit, vLR, fluxes);

  fluxes[1] *= Anorm;
  fluxes[2] *= Anorm;
  fluxes[3] *= Anorm;
  fluxes[4] *= Anorm;
}

/**
 * @brief Update the fluxes for the particle with the given contributions,
 * assuming the particle is to the left of the interparticle interface.
 *
 * @param p Particle.
 * @param fluxes Fluxes accross the interface.
 * @param dx Distance between the particles that share the interface.
 * @param dt Time step for the flux exchange.
 */
__attribute__((always_inline)) INLINE static void hydro_part_update_fluxes_left(
    struct part* restrict p, const float* fluxes, const float* dx,
    const float dt) {

  p->flux.momentum[0] -= fluxes[1] * dt;
  p->flux.momentum[1] -= fluxes[2] * dt;
  p->flux.momentum[2] -= fluxes[3] * dt;
  p->flux.energy -= fluxes[4] * dt;

#ifndef GIZMO_TOTAL_ENERGY
  p->flux.energy += fluxes[1] * p->v[0] * dt;
  p->flux.energy += fluxes[2] * p->v[1] * dt;
  p->flux.energy += fluxes[3] * p->v[2] * dt;
#endif
}

/**
 * @brief Update the fluxes for the particle with the given contributions,
 * assuming the particle is to the right of the interparticle interface.
 *
 * @param p Particle.
 * @param fluxes Fluxes accross the interface.
 * @param dx Distance between the particles that share the interface.
 * @param dt Time step for the flux exchange.
 */
__attribute__((always_inline)) INLINE static void
hydro_part_update_fluxes_right(struct part* restrict p, const float* fluxes,
                               const float* dx, const float dt) {

  p->flux.momentum[0] += fluxes[1] * dt;
  p->flux.momentum[1] += fluxes[2] * dt;
  p->flux.momentum[2] += fluxes[3] * dt;
  p->flux.energy += fluxes[4] * dt;

#ifndef GIZMO_TOTAL_ENERGY
  p->flux.energy -= fluxes[1] * p->v[0] * dt;
  p->flux.energy -= fluxes[2] * p->v[1] * dt;
  p->flux.energy -= fluxes[3] * p->v[2] * dt;
#endif
}

/**
 * @brief Get the drift term for the density based on the mass flux.
 *
 * @param mass_flux Current mass flux for the particle.
 * @param dt Drift time step (in co-moving units).
 * @param volume Volume of the particle.
 * @return 0, since this is Gizmo MFM.
 */
__attribute__((always_inline)) INLINE static float
hydro_gizmo_mfv_density_drift_term(const float mass_flux, const float dt,
                                   const float volume) {
  return 0.0f;
}

/**
 * @brief Add the gravitational contribution to the fluid velocity drift.
 *
 * (MFV only)
 *
 * @param v (drifted) particle velocity.
 * @param fluid_v Fluid velocity.
 * @param v_full (Undrifted) particle velocity.
 * @param a_grav Gravitational acceleration.
 * @param dt_kick_grav Time-step to kick the particle gravitationally.
 */
__attribute__((always_inline)) INLINE static void
hydro_gizmo_mfv_extra_velocity_drift(float* v, const float* restrict fluid_v,
                                     const float* restrict v_full,
                                     const float* restrict a_grav,
                                     float dt_kick_grav) {}

/**
 * @brief Get the term required to update the MFV energy due to the change in
 * gravitational energy.
 *
 * @param dt_kick_corr Time step for the potential energy correction.
 * @param p Particle.
 * @param momentum Momentum of the particle, explicitly requested so that it is
 * clear from the code that the momentum needs to be updated after the call to
 * this function.
 * @param a_grav Gravitational acceleration.
 * @param grav_kick_factor Gravitational kick factor
 * (a_grav * dt + a_grav_mesh * dt_mesh)
 * @return 0, since this is Gizmo MFM.
 */
__attribute__((always_inline)) INLINE static float
hydro_gizmo_mfv_gravity_energy_update_term(const float dt_kick_corr,
                                           const struct part* restrict p,
                                           const float* momentum,
                                           const float* a_grav,
                                           const float* grav_kick_factor) {

  return 0.0f;
}

/**
 * @brief Get the term required to update the MFV mass due to the mass flux.
 *
 * @param mass_flux Mass flux rate.
 * @param dt Time step (in comoving units).
 * @return 0, since this is Gizmo MFM.
 */
__attribute__((always_inline)) INLINE static float
hydro_gizmo_mfv_mass_update_term(const float mass_flux, const float dt) {
  return 0.0f;
}

/**
 * @brief Update the mass of the gpart associated with the given particle after
 * the mass has been updated with the hydrodynamical mass flux.
 *
 * This function does nothing, since this is Gizmo MFM.
 *
 * @param p Particle.
 */
__attribute__((always_inline)) INLINE static void
hydro_gizmo_mfv_update_gpart_mass(struct part* restrict p) {}

#endif /* SWIFT_GIZMO_MFM_HYDRO_FLUX_H */