/* 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; }; static void pdedef(Integer, double, double, const double[], const double[], Integer, const double[], const double[], double[], double[], double[], double[], Integer *, Nag_Comm *); static void bndry1(Integer, Integer, double, const double[], const double[], Integer, const double[], const double[], Integer, double[], Integer *, Nag_Comm *); static void bndry2(Integer, Integer, double, const double[], const double[], Integer, const double[], const double[], Integer, double[], Integer *, Nag_Comm *); static void nmflx1(Integer, double, double, Integer, const double[], const double[], const double[], double[], Integer *, Nag_Comm *, Nag_D03_Save *); static void nmflx2(Integer, double, double, Integer, const double[], const double[], const double[], double[], Integer *, Nag_Comm *, Nag_D03_Save *); static void odedef(Integer, double, Integer, const double[], const double[], Integer, const double[], const double[], const double [], const double[], double[], Integer *, Nag_Comm *); 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), 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) { Vprintf("nag_pde_parab_1d_cd_ode (d03plc) Example Program Results\n"); ex1(); ex2(); return 0; } int ex1(void) { const Integer npde=2, npts=141, ncode=2, nxi=2, neqn=npde*npts+ncode, outpts=8, lrsave=11000, lisave=15700; double tout, ts; Integer exit_status, i, ii, ind, itask, itol, itrace, j, nop; double *algopt=0, *atol=0, *rsave=0, *rtol=0, *u=0, *ue=0, *uout=0, *x=0, *xi=0, *xout=0; Integer *isave=0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; /* 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)) ) { Vprintf("Allocation failure\n"); exit_status = 1; goto END; } Vprintf("\n\nExample 1\n\n"); INIT_FAIL(fail); exit_status = 0; itrace = 0; itol = 1; atol[0] = 1e-5; rtol[0] = 2.5e-4; Vprintf(" npts = %4ld", npts); Vprintf(" atol = %10.3e", atol[0]); Vprintf(" rtol = %10.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) { Vprintf("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; } Vprintf(" t = %6.3f\n\n", ts); Vprintf(" x Approx u1 Exact u1"); Vprintf(" 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) { Vprintf(" %10.4f", xout[i-1]); Vprintf(" %10.4f", UOUT(1,i)); Vprintf(" %10.4f", UE(1,i)); Vprintf(" %10.4f", UOUT(2,i)); Vprintf(" %10.4f\n", UE(2,i)); } 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 (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 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 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)) ) { Vprintf("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 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 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 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 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, neqn=npde*npts+ncode, outpts=8, lisave=neqn+24, lrsave=16392; double d, p, tout, ts, v; Integer exit_status, i, ind, it, itask, itol, itrace, k; double *algopt=0, *atol=0, *rsave=0, *rtol=0, *u=0, *ue=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)) || !(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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } Vprintf("\n\nExample 2\n\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); INIT_FAIL(fail); exit_status = 0; /* 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; Vprintf(" gamma =%6.3f", data.gamma); Vprintf(" elo =%6.3f", data.elo); Vprintf(" ero =%6.3f", data.ero); Vprintf(" rlo =%6.3f", data.rlo); Vprintf(" rro =%6.3f\n\n", data.rro); Vprintf(" npts = %4ld", npts); Vprintf(" atol = %10.3e", atol[0]); Vprintf(" rtol = %10.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; Vprintf(" 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, d03plp, nmflx2, bndry2, 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 nag_pde_parab_1d_cd_ode (d03plc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\n t = %6.3f\n\n", ts); /* Read exact data at output points */ Vscanf(" %*[^\n] "); for (i = 1; i <= 8; ++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 = 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); Vprintf("%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)); } } 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 (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 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 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 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) { Vprintf("Error from nag_pde_parab_1d_euler_osher (d03pvc).\n%s\n", fail.message); } return; }