/* nag_pde_parab_1d_euler_exact (d03pxc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include /* Structure to communicate with user-supplied function arguments */ struct user { double elo, ero, rlo, rro, ulo, uro, gamma; }; static void bndary(Integer, Integer, double, const double[], const double[], Integer, const double[], const double[], Integer, double[], Integer *, Nag_Comm *); static void numflx(Integer, double, double, Integer, const double[], const double[], const double[], double[], Integer *, Nag_Comm *, Nag_D03_Save *); #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, neqn=npde*npts+ncode, lisave=neqn+24, intpts=9, nwkres=npde*(2*npts+3*npde+32)+7*npts+4, lenode=9*neqn+50, mlu=3*npde-1, lrsave=(3*mlu+1)*neqn+nwkres+lenode; double d, p, tout, ts, v; Integer exit_status, i, ind, itask, itol, itrace, k; double *algopt=0, *atol=0, *rtol=0, *u=0, *ue=0, *rsave=0, *x=0, *xi=0; Integer *isave=0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; struct user data; /* 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)) ) { Vprintf("Allocation failure\n"); exit_status = 1; goto END; } INIT_FAIL(fail); exit_status = 0; Vprintf("d03pxc Example Program Results\n"); /* Skip heading in data file */ Vscanf("%*[^\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; d03plc(npde, &ts, tout, d03plp, numflx, bndary, u, npts, x, ncode, d03pek, 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) { Vprintf("Error from d03plc.\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf(" t = %6.3f\n\n", ts); Vprintf(" x APPROX d EXACT d APPROX v EXACT v"); Vprintf(" APPROX p EXACT p\n"); /* Read exact data at output points */ for (i = 1; i <= intpts; ++i) { Vscanf("%lf", &UE(1,i)); Vscanf("%lf", &UE(2,i)); Vscanf("%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); Vprintf(" %8.2e", x[i-1]); Vprintf(" %10.4e", d); Vprintf(" %10.4e", UE(1,k)); Vprintf(" %10.4e", v); Vprintf(" %10.4e", UE(2,k)); Vprintf(" %10.4e", p); Vprintf(" %10.4e\n", UE(3,k)); } Vprintf("\n"); Vprintf(" Number of integration steps in time = %6ld\n", isave[0]); Vprintf(" Number of function evaluations = %6ld\n", isave[1]); Vprintf(" Number of Jacobian evaluations =%6ld\n", isave[2]); Vprintf(" Number of iterations = %6ld\n\n", isave[4]); END: if (algopt) NAG_FREE(algopt); if (atol) NAG_FREE(atol); if (rtol) NAG_FREE(rtol); if (u) NAG_FREE(u); if (ue) NAG_FREE(ue); if (rsave) NAG_FREE(rsave); if (x) NAG_FREE(x); if (xi) NAG_FREE(xi); if (isave) NAG_FREE(isave); return exit_status; } static void 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 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); d03pxc(uleft, uright, data->gamma, tol, niter, flux, saved, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from d03pxc.\n%s\n", fail.message); } return; }