// SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
// SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause

/// @file
/// QFunctions for the `bc_slip` boundary conditions
#include "bc_freestream_type.h"
#include "newtonian_state.h"
#include "newtonian_types.h"
#include "riemann_solver.h"

CEED_QFUNCTION_HELPER int Slip(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
  const NewtonianIdealGasContext newt_ctx = (const NewtonianIdealGasContext)ctx;
  const CeedScalar(*q)[CEED_Q_VLA]        = (const CeedScalar(*)[CEED_Q_VLA])in[0];
  const CeedScalar(*q_data_sur)           = in[2];
  CeedScalar(*v)[CEED_Q_VLA]              = (CeedScalar(*)[CEED_Q_VLA])out[0];
  CeedScalar(*jac_data_sur)               = newt_ctx->is_implicit ? out[1] : NULL;

  const NewtonianIGProperties gas = newt_ctx->gas;

  CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
    const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
    State            s     = StateFromQ(gas, qi, state_var);

    CeedScalar wdetJb, normal[3];
    QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
    wdetJb *= newt_ctx->is_implicit ? -1. : 1.;

    CeedScalar       vel_reflect[3];
    const CeedScalar vel_normal = Dot3(s.Y.velocity, normal);
    for (CeedInt j = 0; j < 3; j++) vel_reflect[j] = s.Y.velocity[j] - 2. * normal[j] * vel_normal;
    const CeedScalar Y_reflect[5] = {s.Y.pressure, vel_reflect[0], vel_reflect[1], vel_reflect[2], s.Y.temperature};
    State            s_reflect    = StateFromY(gas, Y_reflect);

    StateConservative flux = RiemannFlux_HLLC(gas, s, s_reflect, normal);

    CeedScalar Flux[5];
    UnpackState_U(flux, Flux);
    for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];

    if (newt_ctx->is_implicit) StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
  }
  return 0;
}

CEED_QFUNCTION(Slip_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  return Slip(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
}

CEED_QFUNCTION(Slip_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  return Slip(ctx, Q, in, out, STATEVAR_PRIMITIVE);
}

CEED_QFUNCTION(Slip_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  return Slip(ctx, Q, in, out, STATEVAR_ENTROPY);
}

CEED_QFUNCTION_HELPER int Slip_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
  const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
  const CeedScalar(*q_data_sur)     = in[2];
  const CeedScalar(*jac_data_sur)   = in[4];

  CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];

  const NewtonianIdealGasContext newt_ctx = (const NewtonianIdealGasContext)ctx;
  const NewtonianIGProperties    gas      = newt_ctx->gas;

  CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
    CeedScalar wdetJb, normal[3];
    QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
    wdetJb *= newt_ctx->is_implicit ? -1. : 1.;

    CeedScalar qi[5], dqi[5];
    StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi);
    for (int j = 0; j < 5; j++) dqi[j] = dq[j][i];
    State s  = StateFromQ(gas, qi, state_var);
    State ds = StateFromQ_fwd(gas, s, dqi, state_var);

    CeedScalar       vel_reflect[3];
    const CeedScalar vel_normal = Dot3(s.Y.velocity, normal);
    for (CeedInt j = 0; j < 3; j++) vel_reflect[j] = s.Y.velocity[j] - 2. * normal[j] * vel_normal;
    const CeedScalar Y_reflect[5] = {s.Y.pressure, vel_reflect[0], vel_reflect[1], vel_reflect[2], s.Y.temperature};
    State            s_reflect    = StateFromY(gas, Y_reflect);

    CeedScalar       dvel_reflect[3];
    const CeedScalar dvel_normal = Dot3(ds.Y.velocity, normal);
    for (CeedInt j = 0; j < 3; j++) dvel_reflect[j] = ds.Y.velocity[j] - 2. * normal[j] * dvel_normal;
    const CeedScalar dY_reflect[5] = {ds.Y.pressure, dvel_reflect[0], dvel_reflect[1], dvel_reflect[2], ds.Y.temperature};
    State            ds_reflect    = StateFromY_fwd(gas, s_reflect, dY_reflect);

    StateConservative dflux = RiemannFlux_HLLC_fwd(gas, s, ds, s_reflect, ds_reflect, normal);

    CeedScalar dFlux[5];
    UnpackState_U(dflux, dFlux);
    for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j];
  }
  return 0;
}

CEED_QFUNCTION(Slip_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  return Slip_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
}

CEED_QFUNCTION(Slip_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  return Slip_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE);
}

CEED_QFUNCTION(Slip_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
  return Slip_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY);
}