/* nag_pde_parab_1d_cd_ode_remesh (d03psc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.2, 2017.
 */

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

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

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

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

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

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

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

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

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

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

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

  static void NAG_CALL nmflx2(Integer, double, double, Integer,
                              const double[], const double[],
                              const double[], double[], Integer *,
                              Nag_Comm *, Nag_D03_Save *);
#ifdef __cplusplus
}
#endif

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

#define P(I, J)       p[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, K) uout[npde*(intpts*((K) -1)+(J) -1)+(I) -1]

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

  printf("nag_pde_parab_1d_cd_ode_remesh (d03psc) Example Program Results\n");
  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();

  return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0 : 1;
}

int ex1(void)
{
  const Integer iprint = 1;
  const Integer npts = 61, intpts = 7, itype = 1;
  const Integer npde = 1, ncode = 0, nxi = 0, nxfix = 0;
  const Integer neqn = npde * npts + ncode, lisave = 25 + nxfix + neqn;
  const Integer nwkres = npde * (3 * npts + 3 * npde + 32) + 7 * npts + 3;
  const Integer lenode = 11 * neqn + 50, mlu = 3 * npde - 1;
  const Integer lrsave = (3 * mlu + 1) * neqn + nwkres + lenode;
  static double ruser1[5] = { -1.0, -1.0, -1.0, -1.0, -1.0 };
  static Integer iuser1[1] = {1};
  static double xout[7] = { .2, .3, .4, .5, .6, .7, .8 };
  double        con, dxmesh, tout, trmesh, ts, xratio;
  Integer       exit_status = 0;
  Integer       i, ind, ipminf, it, itask, itol, itrace, m, nrmesh;
  Nag_Boolean   remesh;
  double        *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0;
  double        *uout = 0, *x = 0, *xfix = 0, *xi = 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;
  comm.iuser = iuser1;

  /* 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)) ||
      !(uout = NAG_ALLOC(npde * intpts * itype, double)) ||
      !(x = NAG_ALLOC(npts, double)) ||
      !(xfix = NAG_ALLOC(1, double)) ||
      !(xi = NAG_ALLOC(1, double)) ||
      !(isave = NAG_ALLOC(lisave, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = 1;
    goto END;
  }

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

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

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

  /* Initialize mesh */

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

  /* Set remesh parameters */

  remesh = Nag_TRUE;
  nrmesh = 3;
  dxmesh = 0.0;
  trmesh = 0.0;
  con = 2.0 / (npts - 1.0);
  xratio = 1.5;
  ipminf = 0;

  xi[0] = 0.0;
  ind = 0;
  itask = 1;

  for (i = 0; i < 30; ++i)
    algopt[i] = 0.0;

  /* b.d.f. integration */

  algopt[0] = 1.0;
  algopt[12] = 0.005;

  /* Loop over output value of t */

  ts = 0.0;
  for (it = 0; it < 3; ++it) {
    tout = 0.1 * (it + 1);

    /* nag_pde_parab_1d_cd_ode_remesh (d03psc).
     * 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, remeshing, one space variable
     */
    nag_pde_parab_1d_cd_ode_remesh(npde, &ts, tout, pdef1, nmflx1, bndry1,
                                   uvin1, u, npts, x, ncode, NULLFN, nxi, xi,
                                   neqn, rtol, atol, itol, Nag_OneNorm,
                                   Nag_LinAlgBand, algopt, remesh, nxfix,
                                   xfix, nrmesh, dxmesh, trmesh, ipminf,
                                   xratio, con, monit1, rsave, lrsave, isave,
                                   lisave, itask, itrace, 0, &ind,
                                   &comm, &saved, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_pde_parab_1d_cd_ode_remesh (d03psc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

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

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

    /* Interpolate at output points */

    m = 0;
    /* nag_pde_interp_1d_fd (d03pzc). PDEs, spatial interpolation with
     * nag_pde_parab_1d_cd_ode_remesh (d03psc).
     */
    nag_pde_interp_1d_fd(npde, m, u, npts, x, xout, intpts, itype, uout,
                         &fail);

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

    printf(" Approx u ");

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

  if (iprint>0) {
    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(algopt);
  NAG_FREE(atol);
  NAG_FREE(rsave);
  NAG_FREE(rtol);
  NAG_FREE(u);
  NAG_FREE(uout);
  NAG_FREE(x);
  NAG_FREE(xfix);
  NAG_FREE(xi);
  NAG_FREE(isave);

  return exit_status;
}

static void NAG_CALL uvin1(Integer npde, Integer npts, Integer nxi,
                           const double x[], const double xi[], double u[],
                           Integer ncode, double v[], Nag_Comm *comm)
{
  Integer i;

  if (comm->user[0] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback uvin1, first invocation.)\n");
    }
    comm->user[0] = 0.0;
  }
  for (i = 1; i <= npts; ++i) {
    if (x[i - 1] > 0.2 && x[i - 1] <= 0.4) {
      U(1, i) = sin(nag_pi * (5.0 * x[i - 1] - 1.0));
    }
    else {
      U(1, i) = 0.0;
    }
  }
  return;
}

static void NAG_CALL pdef1(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)
{
  if (comm->user[1] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback pdef1, first invocation.)\n");
    }
    comm->user[1] = 0.0;
  }
  P(1, 1) = 1.0;
  c[0] = 0.002;
  d[0] = ux[0];
  s[0] = 0.0;

  return;
}

static void NAG_CALL 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)
{
  /* Zero solution at both boundaries */

  if (comm->user[2] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback bndry1, first invocation.)\n");
    }
    comm->user[2] = 0.0;
  }
  if (ibnd == 0) {
    g[0] = U(1, 1);
  }
  else {
    g[0] = U(1, npts);
  }
  return;
}

static void NAG_CALL monit1(double t, Integer npts, Integer npde,
                            const double x[], const double u[], double fmon[],
                            Nag_Comm *comm)
{
  double h1, h2, h3;
  Integer i;

  if (comm->user[3] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback monit1, first invocation.)\n");
    }
    comm->user[3] = 0.0;
  }
  for (i = 2; i <= npts - 1; ++i) {
    h1 = x[i - 1] - x[i - 2];
    h2 = x[i] - x[i - 1];
    h3 = 0.5 * (x[i] - x[i - 2]);

    /* Second derivatives */

    fmon[i - 1] = fabs(((U(1, i + 1) - U(1, i)) / h2 -
                        (U(1, i) - U(1, i - 1)) / h1) / h3);
  }
  fmon[0] = fmon[1];
  fmon[npts - 1] = fmon[npts - 2];

  return;
}

static void NAG_CALL 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)
{
  if (comm->user[4] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback nmflx1, first invocation.)\n");
    }
    comm->user[4] = 0.0;
  }
  flux[0] = uleft[0];

  return;
}

int ex2(void)
{
  const Integer iprint = 0;
  const Integer npts = 61, intpts = 7, itype = 1;
  const Integer npde = 1, ncode = 0, nxi = 0, nxfix = 0;
  const Integer neqn = npde * npts + ncode, lisave = 25 + nxfix + neqn;
  const Integer nwkres = npde * (3 * npts + 3 * npde + 32) + 7 * npts + 3;
  const Integer lenode = 11 * neqn + 50, mlu = 3 * npde - 1;
  const Integer lrsave = (3 * mlu + 1) * neqn + nwkres + lenode;
  static double ruser2[5] = { -1.0, -1.0, -1.0, -1.0, -1.0 };
  static Integer iuser2[1] = {0};
  double        con, dxmesh, tout, trmesh, ts, xratio, dx, du, pu, xx;
  Integer       exit_status = 0;
  Integer       i, ind, ipminf, it, itask, itol, itrace, nrmesh;
  Nag_Boolean   remesh;
  double        *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0, *ue = 0;
  double        *uout = 0, *x = 0, *xfix = 0, *xi = 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;
  comm.iuser = iuser2;

  /* 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(npts, double)) ||
      !(ue = NAG_ALLOC(npde * intpts, double)) ||
      !(uout = NAG_ALLOC(npde * intpts * itype, double)) ||
      !(x = NAG_ALLOC(npts, double)) ||
      !(xfix = NAG_ALLOC(1, double)) ||
      !(xi = NAG_ALLOC(1, double)) ||
      !(isave = NAG_ALLOC(lisave, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = 1;
    goto END;
  }

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

  itrace = 0;
  itol = 1;
  atol[0] = 5e-4;
  rtol[0] = 5e-2;
  printf(" npts = %4" NAG_IFMT "", npts);
  printf(" atol = %12.3e", atol[0]);
  printf(" rtol = %12.3e\n\n", rtol[0]);

  /* Initialize mesh */

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

  /* Set remesh parameters */

  remesh = Nag_TRUE;
  nrmesh = 5;
  dxmesh = 0.0;
  trmesh = 0.0;
  con = 1.0 / (npts - 1.0);
  xratio = 1.5;
  ipminf = 0;
  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.0025;

  printf(" Track centre point, x,  of moving front\n\n");
  printf("    t        x        Exact x\n\n");
  ts = 0.0;
  tout = 0.0;
  for (it = 0; it < 10; ++it) {
    tout = tout + 0.02;

    /* nag_pde_parab_1d_cd_ode_remesh (d03psc), see above. */
    nag_pde_parab_1d_cd_ode_remesh(npde, &ts, tout, pdef2, nmflx2, bndry2,
                                   uvin2, u, npts, x, ncode, NULLFN, nxi, xi,
                                   neqn, rtol, atol, itol, Nag_OneNorm,
                                   Nag_LinAlgBand, algopt, remesh, nxfix,
                                   xfix, nrmesh, dxmesh, trmesh, ipminf,
                                   xratio, con, monit2, rsave, lrsave, isave,
                                   lisave, itask, itrace, 0, &ind,
                                   &comm, &saved, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_pde_parab_1d_cd_ode_remesh (d03psc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    /* Find point x at which u = 0.5 */
    i = 0;
    while (u[i]>0.5) {
      i++;
    }
    du = u[i-1] - u[i];
    dx = x[i] - x[i-1];
    pu = (1.0 - u[i-1])/du;
    xx = x[i-1] + pu*dx;
    
    /* Check against exact solution */

    printf(" %6.2f", ts);
    printf("  %8.2f", xx);
    printf("   %8.3f\n", ts + 21.0/200.0);
  }

  if (iprint>0) {
    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(algopt);
  NAG_FREE(atol);
  NAG_FREE(rsave);
  NAG_FREE(rtol);
  NAG_FREE(u);
  NAG_FREE(ue);
  NAG_FREE(uout);
  NAG_FREE(x);
  NAG_FREE(xfix);
  NAG_FREE(xi);
  NAG_FREE(isave);

  return exit_status;
}

static void NAG_CALL uvin2(Integer npde, Integer npts, Integer nxi,
                           const double x[], const double xi[], double u[],
                           Integer ncode, double v[], Nag_Comm *comm)
{
  double t;

  if (comm->user[0] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback uvin2, first invocation.)\n");
    }
    comm->user[0] = 0.0;
  }
  t = 0.0;
  exact(t, u, x, npde, npts);

  return;
}

static void NAG_CALL pdef2(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)
{
  if (comm->user[1] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback pdef2, first invocation.)\n");
    }
    comm->user[1] = 0.0;
  }
  P(1, 1) = 1.0;
  c[0] = 0.0;
  d[0] = 0.0;
  s[0] = -100.0 * u[0] * (u[0] - 1.0) * (u[0] - 0.5);

  return;
}

static void NAG_CALL 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)
{
  /* Solution known to be constant at both boundaries */

  double ue[1];

  if (comm->user[2] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback bndry2, first invocation.)\n");
    }
    comm->user[2] = 0.0;
  }
  if (ibnd == 0) {
    exact(t, ue, &x[0], npde, 1);
    g[0] = UE(1, 1) - U(1, 1);
  }
  else {
    exact(t, ue, &x[npts - 1], npde, 1);
    g[0] = UE(1, 1) - U(1, npts);
  }

  return;
}

static void NAG_CALL 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)
{
  if (comm->user[3] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback nmflx2, first invocation.)\n");
    }
    comm->user[3] = 0.0;
  }
  flux[0] = uleft[0];

  return;
}

static void NAG_CALL monit2(double t, Integer npts, Integer npde,
                            const double x[], const double u[],
                            double fmon[], Nag_Comm *comm)
{
  static double xa = 0.0;
  double h1, ux, uxmax, xl, xmax, xr;
  Integer i;

  if (comm->user[4] == -1.0) {
    if (comm->iuser[0] > 0) {
      printf("(User-supplied callback monit2, first invocation.)\n");
    }
    comm->user[4] = 0.0;
  }
  /* Locate shock */

  uxmax = 0.0;
  xmax = 0.0;
  for (i = 2; i <= npts - 1; ++i) {
    h1 = x[i - 1] - x[i - 2];
    ux = fabs((U(1, i) - U(1, i - 1)) / h1);
    if (ux > uxmax) {
      uxmax = ux;
      xmax = x[i - 1];
    }
  }

  xa = 7.0/60.0;
  
  xl = xmax - xa;
  xr = xmax + xa;

  /* Assign monitor function */

  for (i = 0; i < npts; ++i) {
    if (x[i] > xl && x[i] < xr) {
      fmon[i] = 1.0 + cos(nag_pi * (x[i] - xmax) / xa);
    }
    else {
      fmon[i] = 0.0;
    }
  }
  return;
}

static void exact(double t, double *u, const double *x, Integer npde,
                  Integer npts)
{
  /* Exact solution (for comparison and b.c. purposes) */

  double del, psi, rm, rn, s;
  Integer i;
  s = 0.1;
  del = 0.01;
  rm = -1.0 / del;
  rn = s / del + 1.0;

  for (i = 1; i <= npts; ++i) {
    psi = x[i - 1] - t;
    if (psi < s) {
      U(1, i) = 1.0;
    }
    else if (psi > del + s) {
      U(1, i) = 0.0;
    }
    else {
      U(1, i) = rm * psi + rn;
    }
  }
  return;
}