/* nag_pde_parab_1d_euler_hll (d03pwc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include #include /* 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 UE(I, J) ue[npde*((J) -1)+(I) -1] #define U(I, J) u[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, i, ind, itask, itol, itrace, k; double *algopt = 0, *atol = 0, *rtol = 0, *rsave, *u = 0, *ue = 0; double *x = 0, *xi = 0; Integer *isave; NagError fail; Nag_Comm comm; Nag_D03_Save saved; struct user data; INIT_FAIL(fail); exit_status = 0; /* 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; } printf( "nag_pde_parab_1d_euler_hll (d03pwc) Example Program Results\n"); /* 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: 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 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; INIT_FAIL(fail); /* nag_pde_parab_1d_euler_hll (d03pwc). * Modified HLL 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_hll(uleft, uright, data->gamma, flux, saved, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_pde_parab_1d_euler_hll (d03pwc).\n%s\n", fail.message); } return; }