/* nag_pde_parab_1d_cd (d03pfc) 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>

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

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

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

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

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

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


static void exact(double, double *, Integer, const double *, Integer);

int main(void)
{
  Vprintf("d03pfc Example Program Results\n");
  ex1();
  ex2();
  return 0;
}


#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, lisave=npde*npts+24, 
    lrsave=(11+9*npde)*npde*npts+(32+3*npde)*npde+7*npts+54;
  Integer exit_status, 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;
  
  /* 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)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    
  Vprintf("\n\nExample 1\n\n\n");
  INIT_FAIL(fail);
  exit_status = 0;

  itrace = 0;
  acc[0] = 1.0e-4;
  acc[1] = 1.0e-5;
  tsmax = 0.0;
    
  Vprintf(" npts = %4ld acc[0] = %10.3e acc[1] = %10.3e\n\n", 
          npts, acc[0], acc[1]);
  Vprintf("        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;

      d03pfc(npde, &ts, tout, d03pfp, numflx1, bndary1, u, npts, 
             x, acc, tsmax, rsave, lrsave, isave, lisave, itask, 
             itrace, 0, &ind, &comm, &saved, &fail);

      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from 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];
        }

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


static void 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, pi;
  double ue[2];

  pi = nag_pi;

  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 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*(-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] - uright[1] + uleft[1]*3.);

  return;
}

  
static void 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, 
    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;
  
  /* 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)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    
  Vprintf("\n\nExample 2\n\n\n");
  INIT_FAIL(fail);

  itrace = 0;
  acc[0] = 1e-5;
  acc[1] = 1e-5;
    
  Vprintf(" npts = %4ld  acc[0] = %10.3e acc[1] = %10.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];

  Vprintf(" x    ");

  for (i = 0; i < 7; ++i) 
    {
      Vprintf("%9.4f", xout[i]);
      Vprintf((i+1)%7 == 0 || i == 6 ?"\n":"");
    }   
  Vprintf("\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;

      d03pfc(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)
        {
          Vprintf("Error from d03pfc.\n%s\n", fail.message);
          exit_status = 1;
          goto END;
        }

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

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


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


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