/* nag_pde_parab_1d_euler_exact (d03pxc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 7, 2001.
*/
#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, rlo, rro, ulo, uro, gamma;
};
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL bndary(Integer, Integer, double, const double[],
const double[], Integer, const double[],
const double[], Integer, double[], Integer *,
Nag_Comm *);
static void NAG_CALL numflx(Integer, double, double, Integer, const double[],
const double[], const double[], double[],
Integer *, Nag_Comm *, Nag_D03_Save *);
#ifdef __cplusplus
}
#endif
#define U(I, J) u[npde*((J) -1)+(I) -1]
#define UE(I, J) ue[npde*((J) -1)+(I) -1]
int main(void)
{
const Integer npde = 3, npts = 141, ncode = 0, nxi = 0;
const Integer neqn = npde*npts+ncode, lisave = neqn+24, intpts = 9;
const Integer nwkres = npde*(2*npts+3*npde+32)+7*npts+4, lenode = 9*neqn+50;
const Integer mlu = 3*npde-1, lrsave = (3*mlu+1)*neqn+nwkres+lenode;
double d, p, tout, ts, v;
Integer exit_status = 0, i, ind, itask, itol, itrace, k;
double *algopt = 0, *atol = 0, *rtol = 0, *u = 0, *ue = 0, *rsave = 0;
double *x = 0, *xi = 0;
Integer *isave = 0;
NagError fail;
Nag_Comm comm;
Nag_D03_Save saved;
struct user data;
INIT_FAIL(fail);
printf(
"nag_pde_parab_1d_euler_exact (d03pxc) Example Program Results\n");
/* Allocate memory */
if (!(algopt = NAG_ALLOC(30, double)) ||
!(atol = NAG_ALLOC(1, double)) ||
!(rtol = NAG_ALLOC(1, double)) ||
!(u = NAG_ALLOC(npde*npts, double)) ||
!(ue = NAG_ALLOC(npde*intpts, double)) ||
!(rsave = NAG_ALLOC(lrsave, double)) ||
!(x = NAG_ALLOC(npts, double)) ||
!(xi = NAG_ALLOC(1, double)) ||
!(isave = NAG_ALLOC(lisave, Integer)))
{
printf("Allocation failure\n");
exit_status = 1;
goto END;
}
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Problem parameters */
data.gamma = 1.4;
data.rlo = 5.99924;
data.rro = 5.99242;
data.ulo = 5.99924*19.5975;
data.uro = -5.99242*6.19633;
data.elo = 460.894/(data.gamma-1.0) + 0.5*data.rlo*19.5975*19.5975;
data.ero = 46.095 /(data.gamma-1.0) + 0.5*data.rro*6.19633*6.19633;
comm.p = (Pointer)&data;
/* Initialise mesh */
for (i = 0; i < npts; ++i) x[i] = i/(npts-1.0);
/* Initial values */
for (i = 1; i <= npts; ++i)
{
if (x[i-1] < 0.5)
{
U(1, i) = data.rlo;
U(2, i) = data.ulo;
U(3, i) = data.elo;
}
else if (x[i-1] == 0.5)
{
U(1, i) = 0.5*(data.rlo + data.rro);
U(2, i) = 0.5*(data.ulo + data.uro);
U(3, i) = 0.5*(data.elo + data.ero);
}
else
{
U(1, i) = data.rro;
U(2, i) = data.uro;
U(3, i) = data.ero;
}
}
itrace = 0;
itol = 1;
atol[0] = 0.005;
rtol[0] = 5e-4;
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;
tout = 0.035;
/* 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, NULLFN, numflx, bndary, u, npts, x,
ncode, 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;
}
printf("\n t = %6.3f\n\n", ts);
printf("%15s%18s%22s\n", "d", "v", "p");
printf("%3s%10s%9s%9s%9s%11s%11s\n", "x", "Approx", "Exact",
"Approx", "Exact", "Approx", "Exact");
/* Read exact data at output points */
for (i = 1; i <= intpts; ++i)
{
scanf("%lf", &UE(1, i));
scanf("%lf", &UE(2, i));
scanf("%lf", &UE(3, i));
}
/* Calculate density, velocity and pressure */
k = 0;
for (i = 15; i <= 127; i += 14)
{
++k;
d = U(1, i);
v = U(2, i)/d;
p = d*(data.gamma-1.0)*(U(3, i)/d - 0.5*v*v);
printf("%4.1f", x[i-1]);
printf("%9.4f", d);
printf("%9.4f", UE(1, k));
printf("%9.4f", v);
printf("%9.4f", UE(2, k));
printf("%13.4e", p);
printf("%13.4e\n", UE(3, k));
}
printf("\n");
printf(" Number of time steps = %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", isave[4]);
END:
NAG_FREE(algopt);
NAG_FREE(atol);
NAG_FREE(rtol);
NAG_FREE(u);
NAG_FREE(ue);
NAG_FREE(rsave);
NAG_FREE(x);
NAG_FREE(xi);
NAG_FREE(isave);
return exit_status;
}
static void NAG_CALL bndary(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)
{
struct user *data = (struct user *) comm->p;
if (ibnd == 0)
{
g[0] = U(1, 1) - data->rlo;
g[1] = U(2, 1) - data->ulo;
g[2] = U(3, 1) - data->elo;
}
else
{
g[0] = U(1, npts) - data->rro;
g[1] = U(2, npts) - data->uro;
g[2] = U(3,
npts) - data->ero;
}
return;
}
static void NAG_CALL numflx(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)
{
struct user *data = (struct user *) comm->p;
NagError fail;
Integer niter = 0;
double tol = 0.0;
INIT_FAIL(fail);
/* nag_pde_parab_1d_euler_exact (d03pxc).
* Exact 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_exact(uleft, uright, data->gamma, tol, niter, flux,
saved, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_pde_parab_1d_euler_exact (d03pxc).\n%s\n",
fail.message);
}
return;
}