/* nag_pde_parab_1d_cd (d03pfc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 7, 2001.
* Mark 7b revised, 2004.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd03.h>
#include <nagx01.h>
#include <math.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_parab_1d_cd (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;
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;
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 = %4ld 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");
/* 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;
/* nag_pde_parab_1d_cd (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_parab_1d_cd(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_parab_1d_cd (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];
}
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));
}
}
printf("\n");
printf(" Number of integration steps in time = %6ld\n", isave[0]);
printf(" Number of function evaluations = %6ld\n", isave[1]);
printf(" Number of Jacobian evaluations = %6ld\n", isave[2]);
printf(" Number of iterations = %6ld\n\n", isave[4]);
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)
{
if (comm->user[1] == -1.0)
{
printf("(User-supplied callback numflx1, first invocation.)\n");
comm->user[1] = 0.0;
}
flux[0] = 0.5*(-1.0*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] - 1.0*uright[1] + 3.0*uleft[1]);
return;
}
static void NAG_CALL 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;
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;
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 = %4ld acc[0] = %12.3e acc[1] = %12.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];
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_parab_1d_cd (d03pfc), see above. */
nag_pde_parab_1d_cd(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_parab_1d_cd (d03pfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf(" t = %6.3f\n", ts);
printf(" 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));
}
printf(" Number of integration steps in time = %6ld\n",
isave[0]);
printf(" Number of function evaluations = %6ld\n", isave[1]);
printf(" Number of Jacobian evaluations = %6ld\n", isave[2]);
printf(" Number of iterations = %6ld\n\n", isave[4]);
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;
}