Example description
/* nag_pde_dim1_parab_convdiff_dae (d03plc) Example Program.
 *
 * Copyright 2020 Numerical Algorithms Group.
 *
 * Mark 27.1, 2020.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>

/* Structure to communicate with user-supplied function arguments */

struct user {
  double elo, ero, gamma, rlo, rro;
};

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL pdedef(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 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 *);

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

static void init1(double, double *, Integer, double *, Integer, Integer);
static void init2(Integer, Integer, double *, double *, Nag_Comm *);
static void exact(double, double *, Integer, const double *, Integer);
static int ex1(void);
static int ex2(void);

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

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

  printf("nag_pde_dim1_parab_convdiff_dae (d03plc) 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) {
  /* Constants */
  const Integer print_statistics = 0;
  const Integer npde = 2, npts = 201, nv = 2, nxi = 2;
  const Integer neqn = npde * npts + nv;
  static double ruser1[4] = {-1.0, -1.0, -1.0, -1.0};

  /* Scalars */
  double tout, ts, errmax, lerr, lwgt;
  Integer exit_status = 0, i, ind, itask, itol, itrace, j;
  Integer nsteps, nfuncs, njacs, niters, ierrmax;
  Integer nwkres, lenode, lisave, lrsave;

  /* Arrays */
  double *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0, *ue = 0;
  double *x = 0, *xi = 0;
  Integer *isave = 0;

  /* Nag Types */
  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");

  /* Allocate memory */
  lisave = 25 * neqn + 24;
  nwkres =
      npde * (2 * npts + 6 * nxi + 3 * npde + 26) + nxi + nv + 7 * npts + 2;
  lenode = 11 * neqn + 50;
  lrsave = 4 * neqn + 11 * neqn / 2 + 1 + nwkres + lenode;
  lisave = lisave * 4;
  lrsave = lrsave * 4;

  if (!(algopt = NAG_ALLOC(30, double)) || !(atol = NAG_ALLOC(1, double)) ||
      !(rsave = NAG_ALLOC(lrsave, double)) || !(rtol = NAG_ALLOC(1, double)) ||
      !(u = NAG_ALLOC(neqn, double)) ||
      !(ue = NAG_ALLOC(npde * npts, double)) ||
      !(x = NAG_ALLOC(npts, double)) || !(xi = NAG_ALLOC(nxi, double)) ||
      !(isave = NAG_ALLOC(lisave, Integer))) {
    printf("Allocation failure\n");
    exit_status = 1;
    goto END;
  }

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

  printf(" Method parameters:\n");
  printf("  Number of mesh points used = %4" NAG_IFMT "\n", npts);
  printf("  Relative tolerance used    = %12.3e\n", rtol[0]);
  printf("  Absolute tolerance used    = %12.3e\n\n", atol[0]);

  /* Initialize mesh */

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

  /* Set initial values */

  ts = 0.0;
  init1(ts, u, npde, x, npts, nv);

  ind = 0;
  itask = 1;

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

  /* BDF integration */

  algopt[0] = 1.0;

  /* Sparse matrix algebra parameters */

  algopt[28] = 0.1;
  algopt[29] = 1.1;

  tout = 0.5;
  /* nag_pde_dim1_parab_convdiff_dae (d03plc).
   * 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, one space variable
   */
  nag_pde_dim1_parab_convdiff_dae(
      npde, &ts, tout, pdedef, nmflx1, bndry1, u, npts, x, nv, odedef, nxi, xi,
      neqn, rtol, atol, itol, Nag_OneNorm, Nag_LinAlgSparse, algopt, rsave,
      lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_pde_dim1_parab_convdiff_dae (d03plc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Check against exact solution */
  exact(tout, ue, npde, &x[0], npts);
  errmax = 0.0;
  for (i = 1; i < npts; i++) {
    lerr = 0.0;
    for (j = 0; j < npde; j++) {
      lwgt = rtol[0] * fabs(ue[i * npde + j]) + atol[0];
      lerr += fabs(u[i * npde + j] - ue[i * npde + j]) / lwgt;
    }
    lerr = lerr / (double)npde;
    errmax = MAX(errmax, lerr);
  }
  ierrmax = 100 * (Integer)(errmax / 100.0) + 100;
  printf("\n Integration Results:\n");
  printf("  Global error is less than %4" NAG_IFMT ""
         " times the local error tolerance.\n",
         ierrmax);

  /* Print integration statistics (reasonably rounded) */
  if (print_statistics == 1) {
    nsteps = 50 * ((isave[0] + 25) / 50);
    nfuncs = 100 * ((isave[1] + 50) / 100);
    njacs = 20 * ((isave[2] + 10) / 20);
    niters = 100 * ((isave[4] + 50) / 100);
    printf("\n Integration Statistics:\n");
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n", "Number of time steps",
           50, nsteps);
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n",
           "Number of function evaluations", 100, nfuncs);
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n",
           "Number of Jacobian evaluations", 20, njacs);
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n", "Number of iterations",
           100, niters);
  }
END:
  NAG_FREE(algopt);
  NAG_FREE(atol);
  NAG_FREE(rsave);
  NAG_FREE(rtol);
  NAG_FREE(u);
  NAG_FREE(ue);
  NAG_FREE(x);
  NAG_FREE(xi);
  NAG_FREE(isave);

  return exit_status;
}

static void NAG_CALL pdedef(Integer npde, double t, double x, const double u[],
                            const double ux[], Integer nv, const double v[],
                            const double vdot[], double p[], double c[],
                            double d[], double s[], Integer *ires,
                            Nag_Comm *comm) {
  Integer i, j;

  if (comm->user[2] == -1.0) {
    /* printf("(User-supplied callback pdedef, first invocation.)\n"); */
    comm->user[2] = 0.0;
  }
  for (i = 1; i <= npde; ++i) {
    c[i - 1] = 1.0;
    d[i - 1] = 0.0;
    s[i - 1] = 0.0;
    for (j = 1; j <= npde; ++j) {
      if (i == j) {
        P(i, j) = 1.0;
      } else {
        P(i, j) = 0.0;
      }
    }
  }
  return;
}

static void NAG_CALL bndry1(Integer npde, Integer npts, double t,
                            const double x[], const double u[], Integer nv,
                            const double v[], const double vdot[], Integer ibnd,
                            double g[], Integer *ires, Nag_Comm *comm) {
  double dudx;
  double *ue = 0;

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

  /* Allocate memory */

  if (!(ue = NAG_ALLOC(npde, double))) {
    printf("Allocation failure\n");
    goto END;
  }

  if (ibnd == 0) {
    exact(t, ue, npde, &x[0], 1);
    g[0] = U(1, 1) + U(2, 1) - UE(1, 1) - UE(2, 1);
    dudx = (U(1, 2) - U(2, 2) - U(1, 1) + U(2, 1)) / (x[1] - x[0]);
    g[1] = vdot[0] - dudx;
  } else {
    exact(t, ue, npde, &x[npts - 1], 1);
    g[0] = U(1, npts) - U(2, npts) - UE(1, 1) + UE(2, 1);
    dudx = (U(1, npts) + U(2, npts) - U(1, npts - 1) - U(2, npts - 1)) /
           (x[npts - 1] - x[npts - 2]);
    g[1] = vdot[1] + 3.0 * dudx;
  }
END:
  NAG_FREE(ue);

  return;
}

static void NAG_CALL nmflx1(Integer npde, double t, double x, Integer nv,
                            const double v[], 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 nmflx1, first invocation.)\n"); */
    comm->user[1] = 0.0;
  }
  flux[0] = 0.5 * (3.0 * uleft[0] - uright[0] + 3.0 * uleft[1] + uright[1]);
  flux[1] = 0.5 * (3.0 * uleft[0] + uright[0] + 3.0 * uleft[1] - uright[1]);
  return;
}

static void NAG_CALL odedef(Integer npde, double t, Integer nv,
                            const double v[], const double vdot[], Integer nxi,
                            const double xi[], const double ucp[],
                            const double ucpx[], const double ucpt[],
                            double r[], Integer *ires, Nag_Comm *comm) {
  if (comm->user[3] == -1.0) {
    /* printf("(User-supplied callback odedef, first invocation.)\n"); */
    comm->user[3] = 0.0;
  }
  if (*ires == -1) {
    r[0] = 0.0;
    r[1] = 0.0;
  } else {
    r[0] = v[0] - UCP(1, 1) + UCP(2, 1);
    r[1] = v[1] - UCP(1, 2) - UCP(2, 2);
  }
  return;
}

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

  double f, g, x1, x2;
  Integer i;

  for (i = 0; i < npts; ++i) {
    x1 = x[i] - 3.0 * t;
    x2 = x[i] + t;
    f = exp(nag_math_pi * x1) * sin(2.0 * nag_math_pi * x1);
    g = exp(-2.0 * nag_math_pi * x2) * cos(2.0 * nag_math_pi * x2);
    u[npde * i] = f + g;
    u[npde * i + 1] = f - g;
  }
  return;
}

static void init1(double t, double *u, Integer npde, double *x, Integer npts,
                  Integer nv) {
  /* Initial solution */

  double f, g, x1, x2;
  Integer i, neqn;

  neqn = npde * npts + nv;
  for (i = 0; i < npts; ++i) {
    x1 = x[i] - 3.0 * t;
    x2 = x[i] + t;
    f = exp(nag_math_pi * x1) * sin(2.0 * nag_math_pi * x1);
    g = exp(-2.0 * nag_math_pi * x2) * cos(2.0 * nag_math_pi * x2);
    u[npde * i] = f + g;
    u[npde * i + 1] = f - g;
  }
  u[neqn - 2] = u[0] - u[1];
  u[neqn - 1] = u[neqn - 3] + u[neqn - 4];

  return;
}

int ex2(void) {
  const Integer print_statistics = 0;
  const Integer npde = 3, npts = 141, nv = 0, nxi = 0;
  const Integer neqn = npde * npts + nv, lisave = neqn + 24;
  const Integer lrsave = 16392;
  static double ruser[2] = {-1.0, -1.0};
  double d, p, tout, ts, v;
  Integer exit_status = 0, i, ind, it, itask, itol, itrace;
  Integer nsteps, nfuncs, njacs, niters;
  double *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0;
  double *x = 0, *xi = 0;
  Integer *isave = 0;
  NagError fail;
  Nag_Comm comm;
  Nag_D03_Save saved;
  struct user data;

  INIT_FAIL(fail);

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

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

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

  /* Problem parameters */

  data.elo = 2.5;
  data.ero = 0.25;
  data.gamma = 1.4;
  data.rlo = 1.0;
  data.rro = 0.125;
  comm.p = (Pointer)&data;

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

  printf(" Problem parameters and initial conditions:\n");
  printf("  gamma          = %5.3f\n", data.gamma);
  printf("      e(x<0.5,0) = %5.3f", data.elo);
  printf("      e(x>0.5,0) = %5.3f\n", data.ero);
  printf("    rho(x<0.5,0) = %5.3f", data.rlo);
  printf("    rho(x>0.5,0) = %5.3f\n\n", data.rro);
  printf(" Method parameters:\n");
  printf("  Number of mesh points used = %4" NAG_IFMT "\n", npts);
  printf("  Relative tolerance used    = %12.3e\n", rtol[0]);
  printf("  Absolute tolerance used    = %12.3e\n\n", atol[0]);

  /* Initialize mesh */

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

  /* Initial values of variables */

  init2(npde, npts, x, u, &comm);

  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.005;
  ts = 0.0;

  printf(" Solution\n%4s%9s%9s%9s%9s\n", "t", "x", "d", "v", "p");
  for (it = 0; it < 2; ++it) {
    tout = 0.1 * (it + 1);

    /* nag_pde_dim1_parab_convdiff_dae (d03plc), see above. */
    nag_pde_dim1_parab_convdiff_dae(
        npde, &ts, tout, NULLFN, nmflx2, bndry2, u, npts, x, nv, NULLFN, nxi,
        xi, neqn, rtol, atol, itol, Nag_TwoNorm, Nag_LinAlgBand, algopt, rsave,
        lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail);

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

    /* Calculate density, velocity and pressure */

    for (i = 1; i <= npts; i += 14) {
      d = U(1, i);
      v = U(2, i) / d;
      p = d * (data.gamma - 1.0) * (U(3, i) / d - 0.5 * v * v);
      if (i == 1) {
        printf("%6.3f  %7.4f  %7.4f  %7.4f  %7.4f\n", ts, x[i - 1], d, v, p);
      } else {
        printf("%6s  %7.4f  %7.4f  %7.4f  %7.4f\n", "", x[i - 1], d, v, p);
      }
    }
    printf("\n");
  }

  /* Print integration statistics (reasonably rounded) */
  if (print_statistics == 1) {
    nsteps = 50 * ((isave[0] + 25) / 50);
    nfuncs = 50 * ((isave[1] + 25) / 50);
    njacs = isave[2];
    niters = isave[4];
    printf("\n Integration Statistics:\n");
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n", "Number of time steps",
           50, nsteps);
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n",
           "Number of function evaluations", 50, nfuncs);
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n",
           "Number of Jacobian evaluations", 1, njacs);
    printf("  %-30s (nearest %3d) = %6" NAG_IFMT "\n", "Number of iterations",
           1, niters);
  }

END:
  NAG_FREE(algopt);
  NAG_FREE(atol);
  NAG_FREE(rsave);
  NAG_FREE(rtol);
  NAG_FREE(u);
  NAG_FREE(x);
  NAG_FREE(xi);
  NAG_FREE(isave);

  return exit_status;
}

static void init2(Integer npde, Integer npts, double *x, double *u,
                  Nag_Comm *comm) {
  Integer i, j;
  struct user *data = (struct user *)comm->p;

  j = 0;
  for (i = 0; i < npts; ++i) {
    if (x[i] < 0.5) {
      u[j] = data->rlo;
      u[j + 1] = 0.0;
      u[j + 2] = data->elo;
    } else if (x[i] == 0.5) {
      u[j] = 0.5 * (data->rlo + data->rro);
      u[j + 1] = 0.0;
      u[j + 2] = 0.5 * (data->elo + data->ero);
    } else {
      u[j] = data->rro;
      u[j + 1] = 0.0;
      u[j + 2] = data->ero;
    }
    j += 3;
  }
  return;
}

static void NAG_CALL bndry2(Integer npde, Integer npts, double t,
                            const double x[], const double u[], Integer nv,
                            const double v[], const double vdot[], Integer ibnd,
                            double g[], Integer *ires, Nag_Comm *comm) {
  struct user *data = (struct user *)comm->p;

  if (comm->user[0] == -1.0) {
    /* printf("(User-supplied callback bndry2, first invocation.)\n"); */
    comm->user[0] = 0.0;
  }
  if (ibnd == 0) {
    g[0] = U(1, 1) - data->rlo;
    g[1] = U(2, 1);
    g[2] = U(3, 1) - data->elo;
  } else {
    g[0] = U(1, npts) - data->rro;
    g[1] = U(2, npts);
    g[2] = U(3, npts) - data->ero;
  }
  return;
}

static void NAG_CALL nmflx2(Integer npde, double t, double x, Integer nv,
                            const double v[], const double uleft[],
                            const double uright[], double flux[], Integer *ires,
                            Nag_Comm *comm, Nag_D03_Save *saved) {
  char solver;
  NagError fail;
  struct user *data = (struct user *)comm->p;

  if (comm->user[1] == -1.0) {
    /* printf("(User-supplied callback nmflx2, first invocation.)\n"); */
    comm->user[1] = 0.0;
  }

  INIT_FAIL(fail);

  solver = 'R';
  if (solver == 'R') {
    /* ROE SCHEME */

    /* nag_pde_dim1_parab_euler_roe (d03puc).
     * Roe's approximate Riemann solver for Euler equations in
     * conservative form, for use with nag_pde_dim1_parab_convdiff
     * (d03pfc), nag_pde_dim1_parab_convdiff_dae (d03plc) and
     * nag_pde_dim1_parab_convdiff_remesh (d03psc)
     */
    nag_pde_dim1_parab_euler_roe(uleft, uright, data->gamma, flux, saved,
                                 &fail);
  } else {
    /* OSHER SCHEME */

    /* nag_pde_dim1_parab_euler_osher (d03pvc).
     * Osher's approximate Riemann solver for Euler equations in
     * conservative form, for use with nag_pde_dim1_parab_convdiff
     * (d03pfc), nag_pde_dim1_parab_convdiff_dae (d03plc) and
     * nag_pde_dim1_parab_convdiff_remesh (d03psc)
     */
    nag_pde_dim1_parab_euler_osher(uleft, uright, data->gamma,
                                   Nag_OsherPhysical, flux, saved, &fail);
  }

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_pde_dim1_parab_euler_osher (d03pvc).\n%s\n",
           fail.message);
  }

  return;
}