/* nag_pde_parab_1d_cd (d03pfc) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 * Mark 7b revised, 2004.
 */

#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);

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

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

  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 = %4ld 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");

  /* Initialise 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 = %6ld\n", isave[0]);
  printf(" Number of function evaluations = %6ld\n", isave[1]);
  printf(" Number of Jacobian evaluations = %6ld\n", isave[2]);
  printf(" Number of iterations = %6ld\n\n", isave[4]);

 END:
  if (acc) NAG_FREE(acc);
  if (rsave) NAG_FREE(rsave);
  if (u) NAG_FREE(u);
  if (ue) NAG_FREE(ue);
  if (x) NAG_FREE(x);
  if (xout) NAG_FREE(xout);
  if (isave) 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 (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)
{
  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;
  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);

  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 = %4ld  acc[0] = %12.3e acc[1] = %12.3e\n\n",
          npts, acc[0], acc[1]);

  /* Initialise 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 = %6ld\n",
          isave[0]);
  printf(" Number of function evaluations = %6ld\n", isave[1]);
  printf(" Number of Jacobian evaluations = %6ld\n", isave[2]);
  printf(" Number of iterations = %6ld\n\n", isave[4]);

 END:
  if (acc) NAG_FREE(acc);
  if (rsave) NAG_FREE(rsave);
  if (u) NAG_FREE(u);
  if (x) NAG_FREE(x);
  if (xout) NAG_FREE(xout);
  if (isave) 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)
{
  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 (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 (x >= 0.0)
    {
      flux [ 0 ] = x * uleft[0];
    }
  else
    {
      flux [ 0 ] = x * uright[0];
    }
  return;
}