/* nag_pde_parab_1d_cd (d03pfc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

#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;
  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;
  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 */

  for (i = 0; i < npts; ++i)
    x[i] = i / (npts - 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) {
      ++nop;
      xout[nop - 1] = x[i];
    }

    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));
    }
  }
  printf("\n");
  printf(" Number of integration steps in time = %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\n", isave[4]);

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)
{
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback numflx1, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  flux[0] =
         0.5 * (-1.0 * uright[0] + 3.0 * uleft[0] + 0.5 * uright[1] +
                1.5 * uleft[1]);
  flux[1] =
         0.5 * (2.0 * uright[0] + 6.0 * uleft[0] - 1.0 * uright[1] +
                3.0 * uleft[1]);

  return;
}

static void NAG_CALL exact(double t, double *u, Integer npde, const double *x,
                           Integer npts)
{
  double x1, x2, pi;
  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;

    U(1, i) = 0.5 * (exp(x1) + exp(x2))
           + 0.25 * (sin(2.0 * pi * (x2 * x2)) - sin(2.0 * pi * (x1 * x1)))
           + 2.0 * t * t - 2.0 * x[i - 1] * t;

    U(2, i) = exp(x2) - exp(x1)
           + 0.5 * (sin(2.0 * pi * (x2 * x2)) + sin(2.0 * pi * (x1 * x1)))
           + x[i - 1] * x[i - 1] + 5.0 * t * t - 2.0 * x[i - 1] * 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;
  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(" t = %6.3f\n", ts);
    printf(" u    %9.4f%9.4f%9.4f%9.4f%9.4f%9.4f%9.4f\n\n", U(1, 1),
           U(1, 4), U(1, 37), U(1, 76), U(1, 112), U(1, 148), U(1, 151));
  }

  printf(" Number of integration steps in time = %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\n", isave[4]);

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;
}