/* nag_pde_dim1_parab_convdiff_remesh (d03psc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*/
#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 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_dim1_parab_convdiff_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_dim1_parab_convdiff_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_dim1_parab_convdiff_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_dim1_parab_convdiff_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_dim1_parab_fd_interp (d03pzc). PDEs, spatial interpolation with
* nag_pde_dim1_parab_convdiff_remesh (d03psc).
*/
nag_pde_dim1_parab_fd_interp(npde, m, u, npts, x, xout, intpts, itype, uout,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_pde_dim1_parab_fd_interp (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_math_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_dim1_parab_convdiff_remesh (d03psc), see above. */
nag_pde_dim1_parab_convdiff_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_dim1_parab_convdiff_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_math_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;
}