/* nag_pde_parab_1d_euler_hll (d03pwc) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd03.h>
#include <nagx01.h>

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

struct user
{
  double elo, ero, rlo, rro, ulo, uro, gamma;
};
static void bndary(Integer, Integer, double, const double[],
                            const double[], Integer, const double[],
                            const double[], Integer, double[], Integer *,
                            Nag_Comm *);
static void numflx(Integer, double, double, Integer, const double[],
                            const double[], const double[], double[], Integer *,
                            Nag_Comm *, Nag_D03_Save *);

#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, neqn=npde*npts+ncode,
    lisave=neqn+24, intpts=9, nwkres=npde*(2*npts+3*npde+32)+7*npts+4,
    lenode=9*neqn+50, 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, *u=0,
    *ue=0, *x=0, *xi=0;
  Integer *isave;
  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)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }

  Vprintf("nag_pde_parab_1d_euler_hll (d03pwc) Example Program Results\n");

  /* Skip heading in data file */

  Vscanf("%*[^\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;

  /* Initialise 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_parab_1d_cd_ode (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_parab_1d_cd_ode(npde, &ts, tout, d03plp, numflx, bndary, u, npts,
                          x, ncode, d03pek, 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)
    {
      Vprintf("Error from nag_pde_parab_1d_cd_ode (d03plc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  Vprintf(" t = %6.3f\n\n", ts);
  Vprintf("    x      APPROX d   EXACT d    APPROX v   EXACT v");
  Vprintf("    APPROX p   EXACT p\n");

  /* Read exact data at output points */

  for (i = 1; i <= intpts; ++i)
    {
      Vscanf("%lf", &UE(1,i));
      Vscanf("%lf", &UE(2,i));
      Vscanf("%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);
      Vprintf(" %8.2e", x[i-1]);
      Vprintf(" %10.4e", d);
      Vprintf(" %10.4e", UE(1,k));
      Vprintf(" %10.4e", v);
      Vprintf(" %10.4e", UE(2,k));
      Vprintf(" %10.4e", p);
      Vprintf(" %10.4e\n", UE(3,k));
    }

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

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

  return exit_status;
}
static void 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 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_parab_1d_euler_hll (d03pwc).
   * Modified HLL Riemann solver for Euler equations in
   * conservative form, for use with nag_pde_parab_1d_cd
   * (d03pfc), nag_pde_parab_1d_cd_ode (d03plc) and
   * nag_pde_parab_1d_cd_ode_remesh (d03psc)
   */
  nag_pde_parab_1d_euler_hll(uleft, uright, data->gamma, flux, saved, &fail);

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

  return;
}