/* nag_pde_parab_1d_cd_ode (d03plc) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.2, 2017.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd03.h>
#include <nagx01.h>
#include <math.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_parab_1d_cd_ode (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_parab_1d_cd_ode (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_parab_1d_cd_ode(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_parab_1d_cd_ode (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_pi * x1) * sin(2.0 * nag_pi * x1);
g = exp(-2.0 * nag_pi * x2) * cos(2.0 * nag_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_pi * x1) * sin(2.0 * nag_pi * x1);
g = exp(-2.0 * nag_pi * x2) * cos(2.0 * nag_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_parab_1d_cd_ode (d03plc), see above. */
nag_pde_parab_1d_cd_ode(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_parab_1d_cd_ode (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 %3" NAG_IFMT ") = %6" NAG_IFMT "\n",
"Number of time steps", 50, nsteps);
printf(" %-30s (nearest %3" NAG_IFMT ") = %6" NAG_IFMT "\n",
"Number of function evaluations", 50, nfuncs);
printf(" %-30s (nearest %3" NAG_IFMT ") = %6" NAG_IFMT "\n",
"Number of Jacobian evaluations", 1, njacs);
printf(" %-30s (nearest %3" NAG_IFMT ") = %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_parab_1d_euler_roe (d03puc).
* Roe's approximate Riemann solver for Euler equations in
* conservative form, for use with nag_pde_parab_1d_cd
* (d03pfc), nag_pde_parab_1d_cd_ode (d03plc) and
* nag_pde_parab_1d_cd_ode_remesh (d03psc)
*/
nag_pde_parab_1d_euler_roe(uleft, uright, data->gamma, flux, saved,
&fail);
}
else {
/* OSHER SCHEME */
/* nag_pde_parab_1d_euler_osher (d03pvc).
* Osher's approximate Riemann solver for Euler equations in
* conservative form, for use with nag_pde_parab_1d_cd
* (d03pfc), nag_pde_parab_1d_cd_ode (d03plc) and
* nag_pde_parab_1d_cd_ode_remesh (d03psc)
*/
nag_pde_parab_1d_euler_osher(uleft, uright, data->gamma,
Nag_OsherPhysical, flux, saved, &fail);
}
if (fail.code != NE_NOERROR) {
printf("Error from nag_pde_parab_1d_euler_osher (d03pvc).\n%s\n",
fail.message);
}
return;
}