/* nag_pde_parab_1d_cd_ode (d03plc) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 */

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

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

struct user
{
  double elo, ero, gamma, rlo, rro;
  FILE   *fpout;
};

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

static void NAG_CALL bndry1(Integer, Integer, double, const double[],
                            const double[], Integer, const double[],
                            const double[], Integer, double[], Integer *,
                            Nag_Comm *);

static void NAG_CALL bndry2(Integer, Integer, double, const double[],
                            const double[], Integer, const double[],
                            const double[], Integer, double[], Integer *,
                            Nag_Comm *);

static void NAG_CALL nmflx1(Integer, double, double, Integer, const double[],
                            const double[], const double[], double[],
                            Integer *, Nag_Comm *, Nag_D03_Save *);

static void NAG_CALL nmflx2(Integer, double, double, Integer, const double[],
                            const double[], const double[], double[],
                            Integer *, Nag_Comm *, Nag_D03_Save *);

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

static void NAG_CALL init1(double, double *, Integer, double *, Integer,
                           Integer);
static void NAG_CALL init2(Integer, Integer, double *, double *, Nag_Comm *);
static void NAG_CALL exact(double, double *, Integer, const double *, Integer);

static int ex1(FILE *fpout, char *outfile);
static int ex2(FILE *fpin, FILE *fpout, char *outfile);

#define P(I, J)    p[npde*((J) -1)+(I) -1]
#define UCP(I, J)  ucp[npde*((J) -1)+(I) -1]
#define UE(I, J)   ue[npde*((J) -1)+(I) -1]
#define U(I, J)    u[npde*((J) -1)+(I) -1]
#define UOUT(I, J) uout[npde*((J) -1)+(I) -1]

int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;
  char     *outfile = 0;
  Integer  exit_status_ex1 = 0;
  Integer  exit_status_ex2 = 0;
  NagError fail;

  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  (void) nag_example_file_io(argc, argv, "-nag_write", &outfile);

  fprintf(fpout, "nag_pde_parab_1d_cd_ode (d03plc) Example Program Results\n");
  exit_status_ex1 = ex1(fpout, outfile);
  exit_status_ex2 = ex2(fpin, fpout, outfile);

  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0: 1;
}

int ex1(FILE *fpout, char *outfile)
{
  const Integer npde = 2, npts = 141, ncode = 2, nxi = 2;
  const Integer neqn = npde*npts+ncode, outpts = 8, lrsave = 11000;
  const Integer lisave = 15700;
  double        tout, ts;
  Integer       exit_status = 0, i, ii, ind, itask, itol, itrace, j, nop;
  double        *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0, *ue = 0;
  double        *uout = 0, *x = 0, *xi = 0, *xout = 0;
  Integer       *isave = 0;
  NagError      fail;
  Nag_Comm      comm;
  Nag_D03_Save  saved;

  INIT_FAIL(fail);

  fprintf(fpout, "\n\nExample 1\n\n");

  comm.p = (Pointer) fpout;

  /* Allocate memory */

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

  itrace = 0;
  itol = 1;
  atol[0] = 1e-5;
  rtol[0] = 2.5e-4;

  fprintf(fpout, " npts = %4ld", npts);
  fprintf(fpout, " atol = %12.3e", atol[0]);
  fprintf(fpout, " rtol = %12.3e\n\n", rtol[0]);

  /* Initialise mesh */

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

  /* Set initial values */

  ts = 0.0;
  init1(ts, u, npde, x, npts, ncode);

  ind = 0;
  itask = 1;

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

  /* Theta integration */

  algopt[0] = 1.0;

  /* Sparse matrix algebra parameters */

  algopt[28] = 0.1;
  algopt[29] = 1.1;

  tout = 0.5;
  /* 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
   */
  if (outfile) fclose(fpout);
  nag_pde_parab_1d_cd_ode(npde, &ts, tout, pdedef, nmflx1, bndry1, u, npts, x,
                          ncode, odedef, nxi, xi, neqn, rtol, atol, itol,
                          Nag_OneNorm, Nag_LinAlgSparse, algopt, rsave, lrsave,
                          isave, lisave, itask, itrace, outfile, &ind, &comm,
                          &saved, &fail);
  if (outfile && !(fpout = fopen(outfile, "a")))
    {
      exit_status = 2;
      goto END;
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_pde_parab_1d_cd_ode (d03plc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Set output points */

  nop = 0;
  for (i = 0; i < npts; i += 20)
    {
      xout[nop] = x[i];
      ++nop;
    }

  fprintf(fpout, " t = %6.3f\n\n", ts);
  fprintf(fpout, "        x        Approx u1   Exact u1");
  fprintf(fpout, "    Approx u2   Exact u2\n\n");

  for (i = 1; i <= nop; ++i)
    {
      ii = (i-1)*20+1;
      j = npde*(ii - 1);
      UOUT(1, i) = u[j];
      UOUT(2, i) = u[j + 1];
    }

  /* Check against exact solution */

  exact(tout, ue, npde, xout, nop);
  for (i = 1; i <= nop; ++i)
    {
      fprintf(fpout, "  %10.4f", xout[i-1]);
      fprintf(fpout, "  %10.4f", UOUT(1, i));
      fprintf(fpout, "  %10.4f", UE(1, i));
      fprintf(fpout, "  %10.4f", UOUT(2, i));
      fprintf(fpout, "  %10.4f\n", UE(2, i));
    }
  fprintf(fpout, "\n");

  fprintf(fpout, " Number of integration steps in time = %6ld\n", isave[0]);
  fprintf(fpout, " Number of function evaluations = %6ld\n", isave[1]);
  fprintf(fpout, " Number of Jacobian evaluations =%6ld\n", isave[2]);
  fprintf(fpout, " Number of iterations = %6ld\n\n", isave[4]);
 END:
  if (algopt) NAG_FREE(algopt);
  if (atol) NAG_FREE(atol);
  if (rsave) NAG_FREE(rsave);
  if (rtol) NAG_FREE(rtol);
  if (u) NAG_FREE(u);
  if (ue) NAG_FREE(ue);
  if (uout) NAG_FREE(uout);
  if (x) NAG_FREE(x);
  if (xi) NAG_FREE(xi);
  if (xout) NAG_FREE(xout);
  if (isave) NAG_FREE(isave);

  return exit_status;
}


static void NAG_CALL pdedef(Integer npde, double t, double x, const double u[],
                            const double ux[], Integer ncode, const double v[],
                            const double vdot[], double p[], double c[],
                            double d[], double s[], Integer *ires,
                            Nag_Comm *comm)
{
  Integer i, j;

  for (i = 1; i <= npde; ++i)
    {
      c[i-1] = 1.0;
      d[i-1] = 0.0;
      s[i-1] = 0.0;
      for (j = 1; j <= npde; ++j)
        {
          if (i == j)
            {
              P(i, j) = 1.0;
            }
          else
            {
              P(i, j) = 0.0;
            }
        }
    }
  return;
}


static void NAG_CALL bndry1(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)
{
  double dudx;
  double *ue = 0;

  FILE   *fpout = (FILE *) comm->p;

  /* Allocate memory */

  if (!(ue = NAG_ALLOC(npde, double)))
    {
      fprintf(fpout, "Allocation failure\n");
      goto END;
    }

  if (ibnd == 0)
    {
      exact(t, ue, npde, &x[0], 1);
      g[0] = U(1, 1) + U(2, 1) - UE(1, 1) - UE(2, 1);
      dudx = (U(1, 2) - U(2, 2) - U(1, 1) + U(2, 1))/(x[1] - x[0]);
      g[1] = vdot[0] - dudx;
    }
  else
    {
      exact(t, ue, npde, &x[npts-1], 1);
      g[0] = U(1, npts) - U(2, npts) - UE(1, 1) + UE(2, 1);
      dudx = (U(1, npts) + U(2, npts) - U(1, npts-1) - U(2, npts-1))/
             (x[npts-1] - x[npts-2]);
      g[1] = vdot[1] + 3.0*dudx;
    }
 END:
  if (ue) NAG_FREE(ue);

  return;
}


static void NAG_CALL nmflx1(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)
{
  flux[0] = 0.5* (3.0*uleft[0] - uright[0] + 3.0*uleft[1] + uright[1]);
  flux[1] = 0.5* (3.0*uleft[0] + uright[0] + 3.0*uleft[1] - uright[1]);
  return;
}

static void NAG_CALL odedef(Integer npde, double t, Integer ncode,
                            const double v[], const double vdot[], Integer nxi,
                            const double xi[], const double ucp[],
                            const double ucpx[], const double ucpt[],
                            double r[], Integer *ires, Nag_Comm *comm)
{
  if (*ires == -1)
    {
      r[0] = 0.0;
      r[1] = 0.0;
    }
  else
    {
      r[0] = v[0] - UCP(1, 1) + UCP(2, 1);
      r[1] = v[1] - UCP(1, 2) - UCP(2, 2);
    }
  return;
}

static void NAG_CALL exact(double t, double *u, Integer npde, const double *x,
                           Integer npts)
{
  /* Exact solution (for comparison and b.c. purposes) */

  double  f, g;
  Integer i;

  for (i = 1; i <= npts; ++i)
    {
      f = exp(nag_pi* (x[i- 1] - 3.0*t))*sin(2.0*nag_pi*(x[i-1] - 3.0*t));
      g = exp(-2.0*nag_pi* (x[i- 1] + t))*cos(2.0*nag_pi*(x[i-1] + t));
      U(1, i) = f + g;
      U(2, i) = f - g;
    }
  return;
}

static void NAG_CALL init1(double t, double *u, Integer npde, double *x,
                           Integer npts, Integer ncode)
{
  /* Initial solution */

  double  f, g;
  Integer i, j, neqn;

  neqn = npde*npts+ ncode;
  j = 0;
  for (i = 0; i < npts; ++i)
    {
      f = exp(nag_pi* (x[i] - 3.0*t))*sin(2.0*nag_pi* (x[i] - 3.0*t));
      g = exp(-2.0*nag_pi* (x[i] + t))*cos(2.0*nag_pi* (x[i] + t));
      u [j] = f + g;
      u [j+1] = f - g;
      j += 2;
    }
  u[neqn-2] = u[0] - u[1];
  u[ neqn- 1] = u[neqn-3] + u[neqn-4];

  return;
}

int ex2(FILE *fpin, FILE *fpout, char *outfile)
{
  const Integer npde = 3, npts = 141, ncode = 0, nxi = 0;
  const Integer neqn = npde*npts+ncode, outpts = 8, lisave = neqn+24;
  const Integer lrsave = 16392;
  double        d, p, tout, ts, v;
  Integer       exit_status = 0, i, ind, it, itask, itol, itrace, k;
  double        *algopt = 0, *atol = 0, *rsave = 0, *rtol = 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);

  fprintf(fpout, "\n\nExample 2\n\n");

  /* Allocate memory */

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

  /* Skip heading in data file */

  fscanf(fpin, "%*[^\n] ");

  /* Problem parameters */

  data.elo = 2.5;
  data.ero = 0.25;
  data.gamma = 1.4;
  data.rlo = 1.0;
  data.rro = 0.125;
  data.fpout = fpout;
  comm.p = (Pointer)&data;

  itrace = 0;
  itol = 1;
  atol[0] = 0.005;
  rtol[0] = 5e-4;

  fprintf(fpout, " gamma =%6.3f", data.gamma);
  fprintf(fpout, "  elo =%6.3f", data.elo);
  fprintf(fpout, "  ero =%6.3f", data.ero);
  fprintf(fpout, "  rlo =%6.3f", data.rlo);
  fprintf(fpout, "  rro =%6.3f\n\n", data.rro);
  fprintf(fpout, " npts = %4ld", npts);
  fprintf(fpout, " atol = %12.3e", atol[0]);
  fprintf(fpout, " rtol = %12.3e\n\n", rtol[0]);

  /* Initialise mesh */

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

  /* Initial values of variables */

  init2(npde, npts, x, u, &comm);

  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;

  fprintf(fpout,
          "    x    APPROX d EXACT d  APPROX v EXACT v  APPROX p EXACT p\n");
  for (it = 0; it < 2; ++it)
    {
      tout = 0.1* (it+1);

      /* nag_pde_parab_1d_cd_ode (d03plc), see above. */
      if (outfile) fclose(fpout);
      nag_pde_parab_1d_cd_ode(npde, &ts, tout, d03plp, nmflx2, bndry2, u, npts,
                              x, ncode, d03pek, nxi, xi, neqn, rtol, atol,
                              itol, Nag_TwoNorm, Nag_LinAlgBand, algopt, rsave,
                              lrsave, isave, lisave, itask, itrace, outfile,
                              &ind, &comm, &saved, &fail);

      if (outfile && !(fpout = fopen(outfile, "a")))
        {
          exit_status = 2;
          goto END;
        }
      if (fail.code != NE_NOERROR)
        {
          fprintf(fpout,
                  "Error from nag_pde_parab_1d_cd_ode (d03plc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      fprintf(fpout, "\n t = %6.3f\n\n", ts);

      /* Read exact data at output points */

      fscanf(fpin, " %*[^\n] ");
      for (i = 1; i <= 8; ++i)
        {
          fscanf(fpin, "%lf", &UE(1, i));
          fscanf(fpin, "%lf", &UE(2, i));
          fscanf(fpin, "%lf", &UE(3, i));
        }

      /* Calculate density, velocity and pressure */

      k = 0;
      for (i = 29; i <= npts- 14; 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);

          fprintf(fpout, "%7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f\n",
                  x[i-1], d, UE(1, k), v, UE(2, k), p, UE(3, k));
        }
    }

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

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

  return exit_status;
}

static void NAG_CALL init2(Integer npde, Integer npts, double *x, double *u,
                           Nag_Comm *comm)
{
  Integer     i, j;
  struct user *data = (struct user *) comm->p;

  j = 0;
  for (i = 0; i < npts; ++i)
    {
      if (x [i] < 0.5)
        {
          u[j] = data->rlo;
          u[j+1] = 0.0;
          u[j+ 2] = data->elo;
        }
      else
      if (x[i] == 0.5)
        {
          u[j] = 0.5*(data->rlo + data->rro);
          u[ j+ 1] = 0.0;
          u[ j+ 2] = 0.5*(data->elo + data->ero);
        }
      else
        {
          u[ j] = data->rro;
          u [ j+ 1] = 0.0;
          u [ j+ 2] = data->ero;
        }
      j += 3;
    }
  return;
}

static void NAG_CALL bndry2(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);
      g[2] = U(3, 1) - data->elo;
    }
  else
    {
      g[0] = U(1, npts) - data->rro;
      g[1] = U(2, npts);
      g[2] = U(3, npts) - data->ero;
    }
  return;
}

static void NAG_CALL nmflx2(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)
{
  char        solver;
  NagError    fail;
  struct user *data = (struct user *) comm->p;

  INIT_FAIL(fail);

  solver = 'R';
  if (solver == 'R')
    {
      /* ROE SCHEME */

      /* nag_pde_parab_1d_euler_roe (d03puc).
       * Roe's approximate 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_roe(uleft, uright, data->gamma, flux, saved,
                                 &fail);
    }
  else
    {
      /* OSHER SCHEME */

      /* nag_pde_parab_1d_euler_osher (d03pvc).
       * Osher's approximate 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_osher(uleft, uright, data->gamma,
                                   Nag_OsherPhysical, flux, saved, &fail);
    }

  if (fail.code != NE_NOERROR)
    {
      fprintf(data->fpout,
              "Error from nag_pde_parab_1d_euler_osher (d03pvc).\n%s\n",
              fail.message);
    }

  return;
}