/* nag_pde_parab_1d_cd (d03pfc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

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

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

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

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

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

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

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

  static void NAG_CALL exact(double, double *, Integer, const double *,
                             Integer);
#ifdef __cplusplus
}
#endif

int main(void)
{
  Integer exit_status_ex1 = 0;
  Integer exit_status_ex2 = 0;

  printf("nag_pde_parab_1d_cd (d03pfc) Example Program Results\n");
  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();
  return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0 : 1;
}

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

int ex1(void)
{
  double        tout, ts, tsmax, dx;
  const Integer npde = 2, npts = 101, outpts = 7, inter = 20;
  const Integer lisave = npde*npts + 24;
  const Integer lrsave = (11 + 9*npde)*npde*npts + (32 + 3*npde)*npde +
                         7*npts + 54;
  static double ruser1[2] = { -1.0, -1.0 };
  Integer       exit_status = 0, i, ind, it, itask, itrace, j, nop;
  Integer       nsteps, nfuncs, njacs, niters;
  double        *acc = 0, *rsave = 0, *u = 0, *ue = 0, *x = 0, *xout = 0;
  Integer       *isave = 0;
  NagError      fail;
  Nag_Comm      comm;
  Nag_D03_Save  saved;

  INIT_FAIL(fail);

  /* For communication with user-supplied functions: */
  comm.user = ruser1;

  printf("\n\nExample 1\n\n\n");

  /* Allocate memory */

  if (!(acc = NAG_ALLOC(2, double)) ||
      !(rsave = NAG_ALLOC(lrsave, double)) ||
      !(u = NAG_ALLOC(npde * npts, double)) ||
      !(ue = NAG_ALLOC(npde * outpts, double)) ||
      !(x = NAG_ALLOC(npts, double)) ||
      !(xout = NAG_ALLOC(outpts, double)) ||
      !(isave = NAG_ALLOC(lisave, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  itrace = 0;
  acc[0] = 1.0e-4;
  acc[1] = 1.0e-5;
  tsmax = 0.0;

  printf(" npts = %4" NAG_IFMT " acc[0] = %12.3e acc[1] = %12.3e\n\n",
         npts, acc[0], acc[1]);
  printf("         x       Approx u   Exact u   Approx v   Exact v\n");

  /* Initialize mesh */

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

  /* Set initial values */

  ts = 0.0;
  exact(ts, u, npde, x, npts);

  ind = 0;
  itask = 1;

  for (it = 1; it <= 2; ++it) {
    tout = 0.1 * it;

    /* nag_pde_parab_1d_cd (d03pfc).
     * General system of convection-diffusion PDEs with source
     * terms in conservative form, method of lines, upwind
     * scheme using numerical flux function based on Riemann
     * solver, one space variable
     */
    nag_pde_parab_1d_cd(npde, &ts, tout, NULLFN, numflx1, bndary1, u, npts,
                        x, acc, tsmax, rsave, lrsave, isave, lisave, itask,
                        itrace, 0, &ind, &comm, &saved, &fail);

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

    /* Set output points */

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

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

    /* Check against exact solution */

    exact(tout, ue, npde, xout, nop);

    for (i = 1; i <= nop; ++i) {
      j = (i - 1) * inter + 1;
      printf("      %9.4f %9.4f %9.4f %9.4f %9.4f\n",
             xout[i - 1], U(1, j), UE(1, i), U(2, j), UE(2, i));
    }
  }
  nsteps = 5*((isave[0]+2)/5);
  nfuncs = 50*((isave[1]+25)/50);
  njacs = 10*((isave[2]+5)/10);
  niters = 50*((isave[4]+25)/50);
  printf("\n");
  printf(" Number of time steps           (nearest 5)  = %6" NAG_IFMT "\n",
         nsteps);
  printf(" Number of function evaluations (nearest 50) = %6" NAG_IFMT "\n",
         nfuncs);
  printf(" Number of Jacobian evaluations (nearest 10) = %6" NAG_IFMT "\n",
         njacs);
  printf(" Number of iterations           (nearest 50) = %6" NAG_IFMT "\n\n",
         niters);

END:
  NAG_FREE(acc);
  NAG_FREE(rsave);
  NAG_FREE(u);
  NAG_FREE(ue);
  NAG_FREE(x);
  NAG_FREE(xout);
  NAG_FREE(isave);

  return exit_status;
}

void NAG_CALL bndary1(Integer npde, Integer npts, double t, const double x[],
                      const double u[], Integer ibnd, double g[],
                      Integer *ires, Nag_Comm *comm)
{
  double c, exu1, exu2;
  double ue[2];

  if (comm->user[0] == -1.0) {
    /* printf("(User-supplied callback bndary1, first invocation.)\n"); */
    comm->user[0] = 0.0;
  }
  if (ibnd == 0) {
    exact(t, ue, npde, &x[0], 1);
    c = (x[1] - x[0]) / (x[2] - x[1]);
    exu1 = (c + 1.0) * U(1, 2) - c * U(1, 3);
    exu2 = (c + 1.0) * U(2, 2) - c * U(2, 3);
    g[0] = 2.0 * U(1, 1) + U(2, 1) - 2.0 * UE(1, 1) - UE(2, 1);
    g[1] = 2.0 * U(1, 1) - U(2, 1) - 2.0 * exu1 + exu2;
  }
  else {
    exact(t, ue, npde, &x[npts - 1], 1);
    c = (x[npts - 1] - x[npts - 2]) / (x[npts - 2] - x[npts - 3]);
    exu1 = (c + 1.0) * U(1, 2) - c * U(1, 3);
    exu2 = (c + 1.0) * U(2, 2) - c * U(2, 3);
    g[0] = 2.0 * U(1, 1) - U(2, 1) - 2.0 * UE(1, 1) + UE(2, 1);
    g[1] = 2.0 * U(1, 1) + U(2, 1) - 2.0 * exu1 - exu2;
  }

  return;
}

static void NAG_CALL numflx1(Integer npde, double t, double x,
                             const double uleft[], const double uright[],
                             double flux[], Integer *ires, Nag_Comm *comm,
                             Nag_D03_Save *saved)
{
  double  ltmp, rtmp;
  if (comm->user[1] == -1.0) {
    /* printf("(User-supplied callback numflx1, first invocation.)\n"); */
    comm->user[1] = 0.0;
  }
  ltmp = 3.0*uleft[0] + 1.5*uleft[1];
  rtmp = uright[0] - 0.5*uright[1];
  flux[0] = 0.5 * (ltmp-rtmp);
  flux[1] = ltmp + rtmp;
  return;
}

static void NAG_CALL exact(double t, double *u, Integer npde, const double *x,
                           Integer npts)
{
  double  x1, x2, pi, px1, px2;
  Integer i;

  pi = nag_pi;

  /* Exact solution (for comparison and b.c. purposes) */

  for (i = 1; i <= npts; ++i) {
    x1 = x[i-1] + t;
    x2 = x[i-1] - 3.0*t;
    px1 = 0.5*sin(2.0*pi*x1*x1);
    px2 = 0.5*sin(2.0*pi*x2*x2);
    U(1, i) = 0.5*(exp(x1) + exp(x2) + px2 - px1) - t*(x1 + x2);
    U(2, i) = exp(x2) - exp(x1) + px1 + px2 + x1*x2 + 8.0*t*t;
  }
  return;
}

int ex2(void)
{
  double tout, ts, tsmax;
  const Integer npde = 1, npts = 151, outpts = 7, lisave = npde * npts + 24;
  const Integer lrsave = (11 + 9*npde)*npde*npts + (32 + 3*npde)*npde +
                         7*npts + 54;
  static double ruser2[3] = { -1.0, -1.0, -1.0 };
  Integer       exit_status = 0, i, ind, it, itask, itrace;
  Integer       nsteps, nfuncs, njacs, niters;
  double        *acc = 0, *rsave = 0, *u = 0, *x = 0, *xout = 0;
  Integer       *isave = 0;
  NagError      fail;
  Nag_Comm      comm;
  Nag_D03_Save  saved;

  INIT_FAIL(fail);

  /* For communication with user-supplied functions: */
  comm.user = ruser2;

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

  /* Allocate memory */

  if (!(acc = NAG_ALLOC(2, double)) || !(rsave = NAG_ALLOC(lrsave, double))
      || !(u = NAG_ALLOC(npde * npts, double))
      || !(x = NAG_ALLOC(npts, double))
      || !(xout = NAG_ALLOC(outpts, double))
      || !(isave = NAG_ALLOC(lisave, Integer))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  itrace = 0;
  acc[0] = 1e-5;
  acc[1] = 1e-5;

  printf(" npts = %4" NAG_IFMT "  acc[0] = %12.3e acc[1] = %12.3e\n\n",
         npts, acc[0], acc[1]);

  /* Initialize mesh */

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

  /* Set initial values */

  for (i = 1; i <= npts; ++i)
    U(1, i) = x[i - 1] + 4.0;

  ind = 0;
  itask = 1;
  tsmax = 0.02;

  /* Set output points */

  xout[0] = x[0];
  xout[1] = x[3];
  xout[2] = x[36];
  xout[3] = x[75];
  xout[4] = x[111];
  xout[5] = x[147];
  xout[6] = x[150];

  printf(" x          ");

  for (i = 0; i < 7; ++i) {
    printf("%9.4f", xout[i]);
    printf((i + 1) % 7 == 0 || i == 6 ? "\n" : "");
  }
  printf("\n");

  /* Loop over output value of t */

  ts = 0.0;
  tout = 1.0;
  for (it = 0; it < 2; ++it) {
    if (it == 1)
      tout = 10.0;

    /* nag_pde_parab_1d_cd (d03pfc), see above. */
    nag_pde_parab_1d_cd(npde, &ts, tout, pdedef, numflx2, bndary2, u, npts,
                        x, acc, tsmax, rsave, lrsave, isave, lisave, itask,
                        itrace, 0, &ind, &comm, &saved, &fail);

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

    printf(" u(t=%6.3f)%9.4f%9.4f%9.4f%9.4f%9.4f%9.4f%9.4f\n\n", ts, U(1, 1),
           U(1, 4), U(1, 37), U(1, 76), U(1, 112), U(1, 148), U(1, 151));
  }

  nsteps = 5*((isave[0]+2)/5);
  nfuncs = 50*((isave[1]+25)/50);
  njacs = 10*((isave[2]+5)/10);
  niters = 50*((isave[4]+25)/50);
  printf("\n");
  printf(" Number of time steps           (nearest 5)  = %6" NAG_IFMT "\n",
         nsteps);
  printf(" Number of function evaluations (nearest 50) = %6" NAG_IFMT "\n",
         nfuncs);
  printf(" Number of Jacobian evaluations (nearest 10) = %6" NAG_IFMT "\n",
         njacs);
  printf(" Number of iterations           (nearest 50) = %6" NAG_IFMT "\n\n",
         niters);

END:
  NAG_FREE(acc);
  NAG_FREE(rsave);
  NAG_FREE(u);
  NAG_FREE(x);
  NAG_FREE(xout);
  NAG_FREE(isave);

  return exit_status;
}

void NAG_CALL pdedef(Integer npde, double t, double x, const double u[],
                     const double ux[], double p[], double c[], double d[],
                     double s[], Integer *ires, Nag_Comm *comm)
{
  if (comm->user[2] == -1.0) {
    /* printf("(User-supplied callback pdedef, first invocation.)\n"); */
    comm->user[2] = 0.0;
  }
  P(1, 1) = 1.0;
  c[0] = 0.01;
  d[0] = ux[0];
  s[0] = u[0];

  return;
}

void NAG_CALL bndary2(Integer npde, Integer npts, double t, const double x[],
                      const double u[], Integer ibnd, double g[],
                      Integer *ires, Nag_Comm *comm)
{
  if (comm->user[0] == -1.0) {
    /* printf("(User-supplied callback bndary2, first invocation.)\n"); */
    comm->user[0] = 0.0;
  }
  if (ibnd == 0) {
    g[0] = U(1, 1) - 3.0;
  }
  else {
    g[0] = U(1, 1) - 5.0;
  }
  return;
}

static void NAG_CALL numflx2(Integer npde, double t, double x,
                             const double uleft[], const double uright[],
                             double flux[], Integer *ires, Nag_Comm *comm,
                             Nag_D03_Save *saved)
{
  if (comm->user[1] == -1.0) {
    /* printf("(User-supplied callback numflx2, first invocation.)\n"); */
    comm->user[1] = 0.0;
  }
  if (x >= 0.0) {
    flux[0] = x * uleft[0];
  }
  else {
    flux[0] = x * uright[0];
  }
  return;
}