/* nag_pde_dim1_parab_convdiff (d03pfc) Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
*
* Mark 29.2, 2023.
*/
#include <math.h>
#include <nag.h>
#include <stdio.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_dim1_parab_convdiff (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, dx;
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;
Integer nsteps, nfuncs, njacs, niters;
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 */
dx = 1.0 / ((double)(npts - 1));
for (i = 0; i < npts; ++i)
x[i] = ((double)i) * dx;
x[npts - 1] = 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_dim1_parab_convdiff (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_dim1_parab_convdiff(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_dim1_parab_convdiff (d03pfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Set output points */
nop = 0;
for (i = 0; i < 101; i += inter) {
xout[nop] = x[i];
++nop;
}
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));
}
}
nsteps = 5 * ((isave[0] + 2) / 5);
nfuncs = 50 * ((isave[1] + 25) / 50);
njacs = 10 * ((isave[2] + 5) / 10);
niters = 50 * ((isave[4] + 25) / 50);
printf("\n");
printf(" Number of time steps (nearest 5) = %6" NAG_IFMT "\n",
nsteps);
printf(" Number of function evaluations (nearest 50) = %6" NAG_IFMT "\n",
nfuncs);
printf(" Number of Jacobian evaluations (nearest 10) = %6" NAG_IFMT "\n",
njacs);
printf(" Number of iterations (nearest 50) = %6" NAG_IFMT "\n\n",
niters);
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) {
double ltmp, rtmp;
if (comm->user[1] == -1.0) {
/* printf("(User-supplied callback numflx1, first invocation.)\n"); */
comm->user[1] = 0.0;
}
ltmp = 3.0 * uleft[0] + 1.5 * uleft[1];
rtmp = uright[0] - 0.5 * uright[1];
flux[0] = 0.5 * (ltmp - rtmp);
flux[1] = ltmp + rtmp;
return;
}
static void NAG_CALL exact(double t, double *u, Integer npde, const double *x,
Integer npts) {
double x1, x2, pi, px1, px2;
Integer i;
pi = nag_math_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;
px1 = 0.5 * sin(2.0 * pi * x1 * x1);
px2 = 0.5 * sin(2.0 * pi * x2 * x2);
U(1, i) = 0.5 * (exp(x1) + exp(x2) + px2 - px1) - t * (x1 + x2);
U(2, i) = exp(x2) - exp(x1) + px1 + px2 + x1 * x2 + 8.0 * t * 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;
Integer nsteps, nfuncs, njacs, niters;
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_dim1_parab_convdiff (d03pfc), see above. */
nag_pde_dim1_parab_convdiff(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_dim1_parab_convdiff (d03pfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf(" u(t=%6.3f)%9.4f%9.4f%9.4f%9.4f%9.4f%9.4f%9.4f\n\n", ts, U(1, 1),
U(1, 4), U(1, 37), U(1, 76), U(1, 112), U(1, 148), U(1, 151));
}
nsteps = 5 * ((isave[0] + 2) / 5);
nfuncs = 50 * ((isave[1] + 25) / 50);
njacs = 10 * ((isave[2] + 5) / 10);
niters = 50 * ((isave[4] + 25) / 50);
printf("\n");
printf(" Number of time steps (nearest 5) = %6" NAG_IFMT "\n",
nsteps);
printf(" Number of function evaluations (nearest 50) = %6" NAG_IFMT "\n",
nfuncs);
printf(" Number of Jacobian evaluations (nearest 10) = %6" NAG_IFMT "\n",
njacs);
printf(" Number of iterations (nearest 50) = %6" NAG_IFMT "\n\n",
niters);
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;
}