/* nag_pde_dim1_parab_euler_hll (d03pwc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.6, 2022.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>

/* Structure to communicate with user-supplied function arguments */

struct user {
  double elo, ero, rlo, rro, ulo, uro, gamma;
};

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL bndary(Integer, Integer, double, const double[],
                            const double[], Integer, const double[],
                            const double[], Integer, double[], Integer *,
                            Nag_Comm *);

static void NAG_CALL numflx(Integer, double, double, Integer, const double[],
                            const double[], const double[], double[], Integer *,
                            Nag_Comm *, Nag_D03_Save *);
#ifdef __cplusplus
}
#endif

#define UE(I, J) ue[npde * ((J)-1) + (I)-1]
#define U(I, J) u[npde * ((J)-1) + (I)-1]

int main(void) {
  const Integer npde = 3, npts = 141, ncode = 0, nxi = 0;
  const Integer neqn = npde * npts + ncode, lisave = neqn + 24, intpts = 9;
  const Integer nwkres = npde * (2 * npts + 3 * npde + 32) + 7 * npts + 4,
                lenode = 9 * neqn + 50;
  const Integer mlu = 3 * npde - 1,
                lrsave = (3 * mlu + 1) * neqn + nwkres + lenode;
  double d, p, tout, ts, v;
  Integer exit_status, i, ind, itask, itol, itrace, k;
  double *algopt = 0, *atol = 0, *rtol = 0, *rsave = 0, *u = 0, *ue = 0;
  double *x = 0, *xi = 0;
  Integer *isave = 0;
  NagError fail;
  Nag_Comm comm;
  Nag_D03_Save saved;
  struct user data;

  INIT_FAIL(fail);

  exit_status = 0;

  /* Allocate memory */

  if (!(algopt = NAG_ALLOC(30, double)) || !(atol = NAG_ALLOC(1, double)) ||
      !(rtol = NAG_ALLOC(1, double)) || !(u = NAG_ALLOC(npde * npts, double)) ||
      !(ue = NAG_ALLOC(npde * intpts, double)) ||
      !(rsave = NAG_ALLOC(lrsave, double)) || !(x = NAG_ALLOC(npts, double)) ||
      !(xi = NAG_ALLOC(1, double)) || !(isave = NAG_ALLOC(lisave, Integer))) {
    printf("Allocation failure\n");
    exit_status = 1;
    goto END;
  }

  printf("nag_pde_dim1_parab_euler_hll (d03pwc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Problem parameters */

  data.gamma = 1.4;
  data.rlo = 5.99924;
  data.rro = 5.99242;
  data.ulo = 5.99924 * 19.5975;
  data.uro = -5.99242 * 6.19633;
  data.elo = 460.894 / (data.gamma - 1.0) + 0.5 * data.rlo * 19.5975 * 19.5975;
  data.ero = 46.095 / (data.gamma - 1.0) + 0.5 * data.rro * 6.19633 * 6.19633;
  comm.p = (Pointer)&data;

  /* Initialize mesh */

  for (i = 0; i < npts; ++i)
    x[i] = i / (npts - 1.0);

  /* Initial values */

  for (i = 1; i <= npts; ++i) {
    if (x[i - 1] < 0.5) {
      U(1, i) = data.rlo;
      U(2, i) = data.ulo;
      U(3, i) = data.elo;
    } else if (x[i - 1] == 0.5) {
      U(1, i) = 0.5 * (data.rlo + data.rro);
      U(2, i) = 0.5 * (data.ulo + data.uro);
      U(3, i) = 0.5 * (data.elo + data.ero);
    } else {
      U(1, i) = data.rro;
      U(2, i) = data.uro;
      U(3, i) = data.ero;
    }
  }

  itrace = 0;
  itol = 1;
  atol[0] = 0.005;
  rtol[0] = 5e-4;
  xi[0] = 0.0;
  ind = 0;
  itask = 1;

  for (i = 0; i < 30; ++i)
    algopt[i] = 0.0;

  /* Theta integration */

  algopt[0] = 2.0;
  algopt[5] = 2.0;
  algopt[6] = 2.0;

  /* Max. time step */

  algopt[12] = 0.005;

  ts = 0.0;
  tout = 0.035;

  /* nag_pde_dim1_parab_convdiff_dae (d03plc).
   * General system of convection-diffusion PDEs with source
   * terms in conservative form, coupled DAEs, method of
   * lines, upwind scheme using numerical flux function based
   * on Riemann solver, one space variable
   */
  nag_pde_dim1_parab_convdiff_dae(
      npde, &ts, tout, NULLFN, numflx, bndary, u, npts, x, ncode, NULLFN, nxi,
      xi, neqn, rtol, atol, itol, Nag_TwoNorm, Nag_LinAlgBand, algopt, rsave,
      lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail);

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_pde_dim1_parab_convdiff_dae (d03plc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\n t = %6.3f\n\n", ts);
  printf("%15s%18s%22s\n", "d", "v", "p");
  printf("%3s%10s%9s%9s%9s%11s%11s\n", "x", "Approx", "Exact", "Approx",
         "Exact", "Approx", "Exact");

  /* Read exact data at output points */

  for (i = 1; i <= intpts; ++i) {
    scanf("%lf", &UE(1, i));
    scanf("%lf", &UE(2, i));
    scanf("%lf", &UE(3, i));
  }

  /* Calculate density, velocity and pressure */

  k = 0;
  for (i = 15; i <= 127; i += 14) {
    ++k;
    d = U(1, i);
    v = U(2, i) / d;
    p = d * (data.gamma - 1.0) * (U(3, i) / d - 0.5 * v * v);
    printf("%4.1f", x[i - 1]);
    printf("%9.4f", d);
    printf("%9.4f", UE(1, k));
    printf("%9.4f", v);
    printf("%9.4f", UE(2, k));
    printf("%13.4e", p);
    printf("%13.4e\n", UE(3, k));
  }

  printf("\n");
  printf(" Number of time steps           = %6" NAG_IFMT "\n", isave[0]);
  printf(" Number of function evaluations = %6" NAG_IFMT "\n", isave[1]);
  printf(" Number of Jacobian evaluations = %6" NAG_IFMT "\n", isave[2]);
  printf(" Number of iterations           = %6" NAG_IFMT "\n", isave[4]);

END:
  NAG_FREE(algopt);
  NAG_FREE(atol);
  NAG_FREE(rtol);
  NAG_FREE(u);
  NAG_FREE(ue);
  NAG_FREE(rsave);
  NAG_FREE(x);
  NAG_FREE(xi);
  NAG_FREE(isave);

  return exit_status;
}

static void NAG_CALL bndary(Integer npde, Integer npts, double t,
                            const double x[], const double u[], Integer ncode,
                            const double v[], const double vdot[], Integer ibnd,
                            double g[], Integer *ires, Nag_Comm *comm) {
  struct user *data = (struct user *)comm->p;

  if (ibnd == 0) {
    g[0] = U(1, 1) - data->rlo;
    g[1] = U(2, 1) - data->ulo;
    g[2] = U(3, 1) - data->elo;
  } else {
    g[0] = U(1, npts) - data->rro;
    g[1] = U(2, npts) - data->uro;
    g[2] = U(3, npts) - data->ero;
  }
  return;
}

static void NAG_CALL numflx(Integer npde, double t, double x, Integer ncode,
                            const double v[], const double uleft[],
                            const double uright[], double flux[], Integer *ires,
                            Nag_Comm *comm, Nag_D03_Save *saved) {
  struct user *data = (struct user *)comm->p;
  NagError fail;

  INIT_FAIL(fail);
  /* nag_pde_dim1_parab_euler_hll (d03pwc).
   * Modified HLL Riemann solver for Euler equations in
   * conservative form, for use with nag_pde_dim1_parab_convdiff
   * (d03pfc), nag_pde_dim1_parab_convdiff_dae (d03plc) and
   * nag_pde_dim1_parab_convdiff_remesh (d03psc)
   */
  nag_pde_dim1_parab_euler_hll(uleft, uright, data->gamma, flux, saved, &fail);

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_pde_dim1_parab_euler_hll (d03pwc).\n%s\n",
           fail.message);
  }

  return;
}