/* nag_pde_parab_1d_cd (d03pfc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include int ex1(void); int ex2(void); static void pdedef(Integer, double, double, const double[], const double[], double[], double[], double[], double[], Integer *, Nag_Comm *); static void bndary1(Integer, Integer, double, const double[], const double[], Integer, double[], Integer *, Nag_Comm *); static void bndary2(Integer, Integer, double, const double[], const double[], Integer, double[], Integer *, Nag_Comm *); static void numflx1(Integer, double, double, const double[], const double[], double[], Integer *, Nag_Comm *, Nag_D03_Save *); static void numflx2(Integer, double, double, const double[], const double[], double[], Integer *, Nag_Comm *, Nag_D03_Save *); static void exact(double, double *, Integer, const double *, Integer); int main(void) { Vprintf("d03pfc Example Program Results\n"); ex1(); ex2(); return 0; } #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, lisave=npde*npts+24, lrsave=(11+9*npde)*npde*npts+(32+3*npde)*npde+7*npts+54; Integer exit_status, 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; /* 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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } Vprintf("\n\nExample 1\n\n\n"); INIT_FAIL(fail); exit_status = 0; itrace = 0; acc[0] = 1.0e-4; acc[1] = 1.0e-5; tsmax = 0.0; Vprintf(" npts = %4ld acc[0] = %10.3e acc[1] = %10.3e\n\n", npts, acc[0], acc[1]); Vprintf(" 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; d03pfc(npde, &ts, tout, d03pfp, numflx1, bndary1, u, npts, x, acc, tsmax, rsave, lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from 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]; } Vprintf("\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; Vprintf(" %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)); } } 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 (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; } static void 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, pi; double ue[2]; pi = nag_pi; 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 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*(-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] - uright[1] + uleft[1]*3.); return; } static void 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, 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; /* 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)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } Vprintf("\n\nExample 2\n\n\n"); INIT_FAIL(fail); itrace = 0; acc[0] = 1e-5; acc[1] = 1e-5; Vprintf(" npts = %4ld acc[0] = %10.3e acc[1] = %10.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]; Vprintf(" x "); for (i = 0; i < 7; ++i) { Vprintf("%9.4f", xout[i]); Vprintf((i+1)%7 == 0 || i == 6 ?"\n":""); } Vprintf("\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; d03pfc(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) { Vprintf("Error from d03pfc.\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf(" t = %6.3f\n", ts); Vprintf(" 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)); } 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 (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; } static void 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; } static void 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 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; }