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

#include <stdio.h>
#include <nag.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;
};
static void pdedef(Integer, double, double, const double[],
                            const double[], Integer, const double[],
                            const double[], double[], double[], double[],
                            double[], Integer *, Nag_Comm *);
static void bndry1(Integer, Integer, double, const double[],
                            const double[], Integer, const double[],
                            const double[], Integer, double[], Integer *,
                            Nag_Comm *);
static void bndry2(Integer, Integer, double, const double[],
                            const double[], Integer, const double[],
                            const double[], Integer, double[], Integer *,
                            Nag_Comm *);
static void nmflx1(Integer, double, double, Integer, const double[],
                            const double[], const double[], double[], Integer *,
                            Nag_Comm *, Nag_D03_Save *);
static void nmflx2(Integer, double, double, Integer, const double[],
                            const double[], const double[], double[], Integer *,
                            Nag_Comm *, Nag_D03_Save *);
static void odedef(Integer, double, Integer, const double[], const double[],
                            Integer, const double[], const double[],
                            const double [], const double[], double[],
                            Integer *, Nag_Comm *);
static void init1(double, double *, Integer, double *, Integer, Integer);
static void init2(Integer, Integer, double *, double *, Nag_Comm *);
static void exact(double, double *, Integer, const double *, Integer);

static int ex1(void), ex2(void);

#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(void)
{
  Vprintf("nag_pde_parab_1d_cd_ode (d03plc) Example Program Results\n");
  ex1();
  ex2();
  return 0;
}


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

  /* 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)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }

  Vprintf("\n\nExample 1\n\n");
  INIT_FAIL(fail);
  exit_status = 0;

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

  Vprintf(" npts = %4ld", npts);
  Vprintf(" atol = %10.3e", atol[0]);
  Vprintf(" rtol = %10.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
   */
  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, 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;
    }

  /* Set output points */

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

  Vprintf(" t = %6.3f\n\n", ts);
  Vprintf("        x        Approx u1   Exact u1");
  Vprintf("    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)
    {
      Vprintf("  %10.4f", xout[i-1]);
      Vprintf("  %10.4f", UOUT(1,i));
      Vprintf("  %10.4f", UE(1,i));
      Vprintf("  %10.4f", UOUT(2,i));
      Vprintf("  %10.4f\n", UE(2,i));
    }
  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 (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 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 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;

  /* Allocate memory */

  if ( !(ue = NAG_ALLOC(npde, double)) )
    {
      Vprintf("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 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 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 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 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(void)
{
  const Integer npde=3, npts=141, ncode=0, nxi=0, neqn=npde*npts+ncode,
    outpts=8, lisave=neqn+24, lrsave=16392;
  double d, p, tout, ts, v;
  Integer exit_status, i, ind, it, itask, itol, itrace, k;
  double *algopt=0, *atol=0, *rsave=0, *rtol=0, *u=0, *ue=0,
    *x=0, *xi=0;
  Integer *isave=0;
  NagError fail;
  Nag_Comm comm;
  Nag_D03_Save saved;
  struct user data;

  /* 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)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  Vprintf("\n\nExample 2\n\n");

  /* Skip heading in data file */

  Vscanf("%*[^\n] ");

  INIT_FAIL(fail);
  exit_status = 0;

  /* Problem parameters */

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

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

  Vprintf(" gamma =%6.3f", data.gamma);
  Vprintf("  elo =%6.3f", data.elo);
  Vprintf("  ero =%6.3f", data.ero);
  Vprintf("  rlo =%6.3f", data.rlo);
  Vprintf("  rro =%6.3f\n\n", data.rro);
  Vprintf(" npts = %4ld", npts);
  Vprintf(" atol = %10.3e", atol[0]);
  Vprintf(" rtol = %10.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;

  Vprintf("    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. */
      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, 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("\n t = %6.3f\n\n", ts);

      /* Read exact data at output points */

      Vscanf(" %*[^\n] ");
      for (i = 1; i <= 8; ++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 = 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);

          Vprintf("%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));

        }
    }

  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 (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 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 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 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)
    {
      Vprintf("Error from nag_pde_parab_1d_euler_osher (d03pvc).\n%s\n",
              fail.message);
    }

  return;
}