/* nag_pde_parab_1d_cd (d03pfc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. * Mark 7b revised, 2004. */ #include #include #include #include #include #include static int ex1(void); static int ex2(void); static void NAG_CALL pdedef(Integer, double, double, const double[], const double[], double[], double[], double[], double[], Integer *, Nag_Comm *); static void NAG_CALL bndary1(Integer, Integer, double, const double[], const double[], Integer, double[], Integer *, Nag_Comm *); static void NAG_CALL bndary2(Integer, Integer, double, const double[], const double[], Integer, double[], Integer *, Nag_Comm *); static void NAG_CALL numflx1(Integer, double, double, const double[], const double[], double[], Integer *, Nag_Comm *, Nag_D03_Save *); static void NAG_CALL numflx2(Integer, double, double, const double[], const double[], double[], Integer *, Nag_Comm *, Nag_D03_Save *); static void NAG_CALL exact(double, double *, Integer, const double *, Integer); int main(void) { Integer exit_status_ex1 = 0; Integer exit_status_ex2 = 0; printf("nag_pde_parab_1d_cd (d03pfc) Example Program Results\n"); exit_status_ex1 = ex1(); exit_status_ex2 = ex2(); return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0 : 1; } #define U(I, J) u[npde*((J) -1)+(I) -1] #define P(I, J) p[npde*((J) -1)+(I) -1] #define UE(I, J) ue[npde*((J) -1)+(I) -1] int ex1(void) { double tout, ts, tsmax; const Integer npde = 2, npts = 101, outpts = 7, inter = 20; const Integer lisave = npde*npts+24; const Integer lrsave = (11+9*npde)*npde*npts+(32+3*npde)*npde+7*npts+54; Integer exit_status = 0, i, ind, it, itask, itrace, j, nop; double *acc = 0, *rsave = 0, *u = 0, *ue = 0, *x = 0, *xout = 0; Integer *isave = 0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; INIT_FAIL(fail); printf("\n\nExample 1\n\n\n"); /* Allocate memory */ if (!(acc = NAG_ALLOC(2, double)) || !(rsave = NAG_ALLOC(lrsave, double)) || !(u = NAG_ALLOC(npde*npts, double)) || !(ue = NAG_ALLOC(npde*outpts, double)) || !(x = NAG_ALLOC(npts, double)) || !(xout = NAG_ALLOC(outpts, double)) || !(isave = NAG_ALLOC(lisave, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } itrace = 0; acc[0] = 1.0e-4; acc[1] = 1.0e-5; tsmax = 0.0; printf(" npts = %4ld acc[0] = %12.3e acc[1] = %12.3e\n\n", npts, acc[0], acc[1]); printf( " x Approx u Exact u Approx v Exact v\n"); /* Initialise mesh */ for (i = 0; i < npts; ++i) x[i] = i/(npts-1.0); /* Set initial values */ ts = 0.0; exact(ts, u, npde, x, npts); ind = 0; itask = 1; for (it = 1; it <= 2; ++it) { tout = 0.1*it; /* nag_pde_parab_1d_cd (d03pfc). * General system of convection-diffusion PDEs with source * terms in conservative form, method of lines, upwind * scheme using numerical flux function based on Riemann * solver, one space variable */ nag_pde_parab_1d_cd(npde, &ts, tout, NULLFN, numflx1, bndary1, u, npts, x, acc, tsmax, rsave, lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_pde_parab_1d_cd (d03pfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Set output points */ nop = 0; for (i = 0; i < 101; i += inter) { ++nop; xout[nop - 1] = x[i]; } printf("\n t = %6.3f\n\n", ts); /* Check against exact solution */ exact(tout, ue, npde, xout, nop); for (i = 1; i <= nop; ++i) { j = (i-1)*inter+1; printf(" %9.4f %9.4f %9.4f %9.4f %9.4f\n", xout[i-1], U(1, j), UE(1, i), U(2, j), 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 (acc) NAG_FREE(acc); if (rsave) NAG_FREE(rsave); if (u) NAG_FREE(u); if (ue) NAG_FREE(ue); if (x) NAG_FREE(x); if (xout) NAG_FREE(xout); if (isave) NAG_FREE(isave); return exit_status; } void NAG_CALL bndary1(Integer npde, Integer npts, double t, const double x[], const double u[], Integer ibnd, double g[], Integer *ires, Nag_Comm *comm) { double c, exu1, exu2; double ue[2]; if (ibnd == 0) { exact(t, ue, npde, &x[0], 1); c = (x[1] - x[0])/(x[2] - x[1]); exu1 = (c + 1.0)*U(1, 2) - c *U(1, 3); exu2 = (c + 1.0)*U(2, 2) - c *U(2, 3); g[0] = 2.0*U(1, 1) + U(2, 1) - 2.0*UE(1, 1) - UE(2, 1); g[1] = 2.0*U(1, 1) - U(2, 1) - 2.0*exu1 + exu2; } else { exact(t, ue, npde, &x[npts-1], 1); c = (x[npts-1] - x[npts - 2])/(x[npts - 2] - x[npts - 3]); exu1 = (c + 1.0)*U(1, 2) - c *U(1, 3); exu2 = (c + 1.0)*U(2, 2) - c *U(2, 3); g[0] = 2.0*U(1, 1) - U(2, 1) - 2.0*UE(1, 1) + UE(2, 1); g[1] = 2.0*U(1, 1) + U(2, 1) - 2.0*exu1 - exu2; } return; } static void NAG_CALL numflx1(Integer npde, double t, double x, const double uleft[], const double uright[], double flux[], Integer *ires, Nag_Comm *comm, Nag_D03_Save *saved) { flux[0] = 0.5*(-1.0*uright[0] + 3.0*uleft[0] + 0.5*uright[1] + 1.5*uleft[1]); flux[1] = 0.5*(2.0*uright[0] + 6.0*uleft[0] - 1.0*uright[1] + 3.0*uleft[1]); return; } static void NAG_CALL exact(double t, double *u, Integer npde, const double *x, Integer npts) { double x1, x2, pi; Integer i; pi = nag_pi; /* Exact solution (for comparison and b.c. purposes) */ for (i = 1; i <= npts; ++i) { x1 = x[i-1] + t; x2 = x[i-1] - 3.0*t; U(1, i) = 0.5*(exp(x1) + exp(x2)) + 0.25*(sin(2.0*pi*(x2*x2)) - sin(2.0*pi*(x1*x1))) + 2.0*t*t - 2.0*x[i-1]*t; U(2, i) = exp(x2) - exp(x1) + 0.5*(sin(2.0*pi*(x2*x2)) + sin(2.0*pi*(x1*x1))) + x[i-1]* x[i-1] + 5.0*t*t - 2.0*x[i-1]*t; } return; } int ex2(void) { double tout, ts, tsmax; const Integer npde = 1, npts = 151, outpts = 7, lisave = npde*npts+24; const Integer lrsave = (11+9*npde)*npde*npts+(32+3*npde)*npde+7*npts+54; Integer exit_status = 0, i, ind, it, itask, itrace; double *acc = 0, *rsave = 0, *u = 0, *x = 0, *xout = 0; Integer *isave = 0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; INIT_FAIL(fail); printf("\n\nExample 2\n\n\n"); /* Allocate memory */ if (!(acc = NAG_ALLOC(2, double)) || !(rsave = NAG_ALLOC(lrsave, double)) || !(u = NAG_ALLOC(npde*npts, double)) || !(x = NAG_ALLOC(npts, double)) || !(xout = NAG_ALLOC(outpts, double)) || !(isave = NAG_ALLOC(lisave, Integer)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } itrace = 0; acc[0] = 1e-5; acc[1] = 1e-5; printf(" npts = %4ld acc[0] = %12.3e acc[1] = %12.3e\n\n", npts, acc[0], acc[1]); /* Initialise mesh */ for (i = 0; i < npts; ++i) x[i] = -1.0 + 2.0*i/ (npts-1.0); /* Set initial values */ for (i = 1; i <= npts; ++i) U(1, i) = x[i-1] + 4.0; ind = 0; itask = 1; tsmax = 0.02; /* Set output points */ xout[0] = x[0]; xout[1] = x[3]; xout[2] = x[36]; xout[3] = x[75]; xout[4] = x[111]; xout[5] = x[147]; xout[6] = x[150]; printf(" x "); for (i = 0; i < 7; ++i) { printf("%9.4f", xout[i]); printf((i+1)%7 == 0 || i == 6?"\n":""); } printf("\n"); /* Loop over output value of t */ ts = 0.0; tout = 1.0; for (it = 0; it < 2; ++it) { if (it == 1) tout = 10.0; /* nag_pde_parab_1d_cd (d03pfc), see above. */ nag_pde_parab_1d_cd(npde, &ts, tout, pdedef, numflx2, bndary2, u, npts, x, acc, tsmax, rsave, lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_pde_parab_1d_cd (d03pfc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" t = %6.3f\n", ts); printf(" u %9.4f%9.4f%9.4f%9.4f%9.4f%9.4f%9.4f\n\n", U(1, 1), U(1, 4), U(1, 37), U(1, 76), U(1, 112), U(1, 148), U(1, 151)); } 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 (acc) NAG_FREE(acc); if (rsave) NAG_FREE(rsave); if (u) NAG_FREE(u); if (x) NAG_FREE(x); if (xout) NAG_FREE(xout); if (isave) NAG_FREE(isave); return exit_status; } void NAG_CALL pdedef(Integer npde, double t, double x, const double u[], const double ux[], double p[], double c[], double d[], double s[], Integer *ires, Nag_Comm *comm) { P(1, 1) = 1.0; c[0] = 0.01; d[0] = ux[0]; s[0] = u[0]; return; } void NAG_CALL bndary2(Integer npde, Integer npts, double t, const double x[], const double u[], Integer ibnd, double g[], Integer *ires, Nag_Comm *comm) { if (ibnd == 0) { g [ 0] = U(1, 1) - 3.0; } else { g [ 0] = U(1, 1) - 5.0; } return; } static void NAG_CALL numflx2(Integer npde, double t, double x, const double uleft [], const double uright [], double flux[], Integer *ires, Nag_Comm *comm, Nag_D03_Save *saved) { if (x >= 0.0) { flux [ 0 ] = x * uleft[0]; } else { flux [ 0 ] = x * uright[0]; } return; }