/* nag_pde_parab_1d_cd_ode (d03plc) 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, 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 NAG_CALL init1(double, double *, Integer, double *, Integer, Integer); static void NAG_CALL init2(Integer, Integer, double *, double *, Nag_Comm *); static void NAG_CALL 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 UE(I, J) ue[npde*((J) -1)+(I) -1] #define U(I, J) u[npde*((J) -1)+(I) -1] #define UOUT(I, J) uout[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) { const Integer npde = 2, npts = 141, ncode = 2, nxi = 2; const Integer neqn = npde*npts+ncode, outpts = 8, lrsave = 11000; const Integer lisave = 15700; double tout, ts; Integer exit_status = 0, i, ii, ind, itask, itol, itrace, j, nop; double *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0, *ue = 0; double *uout = 0, *x = 0, *xi = 0, *xout = 0; Integer *isave = 0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; INIT_FAIL(fail); printf("\n\nExample 1\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(neqn, double)) || !(ue = NAG_ALLOC(npde*outpts, double)) || !(uout = NAG_ALLOC(npde*outpts, double)) || !(x = NAG_ALLOC(npts, double)) || !(xi = NAG_ALLOC(nxi, double)) || !(xout = NAG_ALLOC(outpts, 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(" npts = %4ld", npts); printf(" atol = %12.3e", atol[0]); printf(" rtol = %12.3e\n\n", rtol[0]); /* Initialise 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, ncode); ind = 0; itask = 1; for (i = 0; i < 30; ++i) algopt[i] = 0.0; /* Theta 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, ncode, 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; } /* Set output points */ nop = 0; for (i = 0; i < npts; i += 20) { xout[nop] = x[i]; ++nop; } printf(" t = %6.3f\n\n", ts); printf(" x Approx u1 Exact u1"); printf(" Approx u2 Exact u2\n\n"); for (i = 1; i <= nop; ++i) { ii = (i-1)*20+1; j = npde*(ii - 1); UOUT(1, i) = u[j]; UOUT(2, i) = u[j + 1]; } /* Check against exact solution */ exact(tout, ue, npde, xout, nop); for (i = 1; i <= nop; ++i) { printf(" %10.4f", xout[i-1]); printf(" %10.4f", UOUT(1, i)); printf(" %10.4f", UE(1, i)); printf(" %10.4f", UOUT(2, i)); printf(" %10.4f\n", 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: if (algopt) NAG_FREE(algopt); if (atol) NAG_FREE(atol); if (rsave) NAG_FREE(rsave); if (rtol) NAG_FREE(rtol); if (u) NAG_FREE(u); if (ue) NAG_FREE(ue); if (uout) NAG_FREE(uout); if (x) NAG_FREE(x); if (xi) NAG_FREE(xi); if (xout) NAG_FREE(xout); if (isave) NAG_FREE(isave); return exit_status; } static void NAG_CALL pdedef(Integer npde, double t, double x, const double u[], const double ux[], Integer ncode, const double v[], const double vdot[], double p[], double c[], double d[], double s[], Integer *ires, Nag_Comm *comm) { Integer i, j; 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 ncode, const double v[], const double vdot[], Integer ibnd, double g[], Integer *ires, Nag_Comm *comm) { double dudx; double *ue = 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: if (ue) NAG_FREE(ue); return; } static void NAG_CALL nmflx1(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) { 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 ncode, 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 (*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 NAG_CALL exact(double t, double *u, Integer npde, const double *x, Integer npts) { /* Exact solution (for comparison and b.c. purposes) */ double f, g; Integer i; for (i = 1; i <= npts; ++i) { f = exp(nag_pi* (x[i- 1] - 3.0*t))*sin(2.0*nag_pi*(x[i-1] - 3.0*t)); g = exp(-2.0*nag_pi* (x[i- 1] + t))*cos(2.0*nag_pi*(x[i-1] + t)); U(1, i) = f + g; U(2, i) = f - g; } return; } static void NAG_CALL init1(double t, double *u, Integer npde, double *x, Integer npts, Integer ncode) { /* Initial solution */ double f, g; Integer i, j, neqn; neqn = npde*npts+ ncode; j = 0; for (i = 0; i < npts; ++i) { f = exp(nag_pi* (x[i] - 3.0*t))*sin(2.0*nag_pi* (x[i] - 3.0*t)); g = exp(-2.0*nag_pi* (x[i] + t))*cos(2.0*nag_pi* (x[i] + t)); u [j] = f + g; u [j+1] = f - g; j += 2; } u[neqn-2] = u[0] - u[1]; u[ neqn- 1] = u[neqn-3] + u[neqn-4]; return; } int ex2(void) { const Integer npde = 3, npts = 141, ncode = 0, nxi = 0; const Integer neqn = npde*npts+ncode, outpts = 8, lisave = neqn+24; const Integer lrsave = 16392; double d, p, tout, ts, v; Integer exit_status = 0, i, ind, it, itask, itol, itrace, k; double *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0, *ue = 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("\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)) || !(ue = NAG_ALLOC(npde*outpts, 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; } /* Skip heading in data file */ scanf("%*[^\n] "); /* 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(" gamma =%6.3f", data.gamma); printf(" elo =%6.3f", data.elo); printf(" ero =%6.3f", data.ero); printf(" rlo =%6.3f", data.rlo); printf(" rro =%6.3f\n\n", data.rro); printf(" npts = %4ld", npts); printf(" atol = %12.3e", atol[0]); printf(" rtol = %12.3e\n\n", rtol[0]); /* Initialise 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( " x APPROX d EXACT d APPROX v EXACT v APPROX p EXACT p\n"); 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, 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); /* Read exact data at output points */ scanf(" %*[^\n] "); for (i = 1; i <= 8; ++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 = 29; i <= npts- 14; 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("%7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f\n", x[i-1], d, UE(1, k), v, UE(2, k), p, UE(3, k)); } } 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: if (algopt) NAG_FREE(algopt); if (atol) NAG_FREE(atol); if (rsave) NAG_FREE(rsave); if (rtol) NAG_FREE(rtol); if (u) NAG_FREE(u); if (ue) NAG_FREE(ue); if (x) NAG_FREE(x); if (xi) NAG_FREE(xi); if (isave) NAG_FREE(isave); return exit_status; } static void NAG_CALL 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 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); 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 ncode, 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; 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; }