/* nag_pde_parab_1d_cd_ode_remesh (d03psc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include static int ex1(void); static int ex2(void); #ifdef __cplusplus extern "C" { #endif static void NAG_CALL uvin1(Integer, Integer, Integer, const double[], const double[], double[], Integer, double[], Nag_Comm *); static void NAG_CALL uvin2(Integer, Integer, Integer, const double[], const double[], double[], Integer, double[], Nag_Comm *); static void NAG_CALL pdef1(Integer, double, double, const double[], const double[], Integer, const double[], const double[], double[], double[], double[], double[], Integer *, Nag_Comm *); static void NAG_CALL pdef2(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 monit1(double, Integer, Integer, const double[], const double[], double[], Nag_Comm *); static void NAG_CALL monit2(double, Integer, Integer, const double[], const double[], double[], 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 *); #ifdef __cplusplus } #endif static void NAG_CALL exact(double, double *, const double *, Integer, Integer); #define P(I, J) p[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, K) uout[npde*(intpts*((K) -1)+(J) -1)+(I) -1] int main(void) { Integer exit_status_ex1 = 0; Integer exit_status_ex2 = 0; printf( "nag_pde_parab_1d_cd_ode_remesh (d03psc) 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 = 1, npts = 61, ncode = 0, nxi = 0, nxfix = 0, itype = 1; const Integer neqn = npde*npts+ncode, intpts = 7, lisave = 25+nxfix+neqn; const Integer nwkres = npde*(3*npts+3*npde+32)+7*npts+3, lenode = 11*neqn+50; const Integer mlu = 3*npde-1, lrsave = (3*mlu+1)*neqn+nwkres+lenode; static double xout[7] = { .2, .3, .4, .5, .6, .7, .8 }; double con, dxmesh, tout, trmesh, ts, xratio; Integer exit_status, i, ind, ipminf, it, itask, itol, itrace, m, nrmesh; Nag_Boolean remesh; double *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0; double *uout = 0, *x = 0, *xfix = 0, *xi = 0; Integer *isave = 0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; INIT_FAIL(fail); exit_status = 0; /* 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)) || !(uout = NAG_ALLOC(npde*intpts*itype, double)) || !(x = NAG_ALLOC(npts, double)) || !(xfix = NAG_ALLOC(1, double)) || !(xi = NAG_ALLOC(1, double)) || !(isave = NAG_ALLOC(lisave, Integer))) { printf("Allocation failure\n"); exit_status = 1; goto END; } printf("\n\nExample 1\n\n"); itrace = 0; itol = 1; atol[0] = 1.0e-4; rtol[0] = 1.0e-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); xfix[0] = 0.0; /* Set remesh parameters */ remesh = Nag_TRUE; nrmesh = 3; dxmesh = 0.0; trmesh = 0.0; con = 2.0/(npts-1.0); xratio = 1.5; ipminf = 0; xi[0] = 0.0; ind = 0; itask = 1; for (i = 0; i < 30; ++i) algopt[i] = 0.0; /* b.d.f. integration */ algopt[0] = 1.0; algopt[12] = 0.005; /* Loop over output value of t */ ts = 0.0; tout = 0.0; for (it = 0; it < 3; ++it) { tout = 0.1*(it+1); /* nag_pde_parab_1d_cd_ode_remesh (d03psc). * 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, remeshing, one space variable */ nag_pde_parab_1d_cd_ode_remesh(npde, &ts, tout, pdef1, nmflx1, bndry1, uvin1, u, npts, x, ncode, NULLFN, nxi, xi, neqn, rtol, atol, itol, Nag_OneNorm, Nag_LinAlgBand, algopt, remesh, nxfix, xfix, nrmesh, dxmesh, trmesh, ipminf, xratio, con, monit1, 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_remesh (d03psc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" t = %6.3f\n", ts); printf(" x "); for (i = 1; i <= intpts; ++i) { printf("%9.4f", xout[i-1]); printf(i%7 == 0 || i == 7?"\n":""); } /* Interpolate at output points */ m = 0; /* nag_pde_interp_1d_fd (d03pzc). PDEs, spatial interpolation with * nag_pde_parab_1d_cd_ode_remesh (d03psc). */ nag_pde_interp_1d_fd(npde, m, u, npts, x, xout, intpts, itype, uout, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_pde_interp_1d_fd (d03pzc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" Approx u "); for (i = 1; i <= intpts; ++i) { printf("%9.4f", UOUT(1, i, 1)); printf(i%7 == 0 || i == 7?"\n":""); } 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 (uout) NAG_FREE(uout); if (x) NAG_FREE(x); if (xfix) NAG_FREE(xfix); if (xi) NAG_FREE(xi); if (isave) NAG_FREE(isave); return exit_status; } static void NAG_CALL uvin1(Integer npde, Integer npts, Integer nxi, const double x[], const double xi[], double u[], Integer ncode, double v[], Nag_Comm *comm) { Integer i; for (i = 1; i <= npts; ++i) { if (x[i-1] > 0.2 && x[i-1] <= 0.4) { U(1, i) = sin(nag_pi*(5.0*x[i-1]-1.0)); } else { U(1, i) = 0.0; } } return; } static void NAG_CALL pdef1(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) { P(1, 1) = 1.0; c[0] = 0.002; d[0] = ux[0]; s[0] = 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) { /* Zero solution at both boundaries */ if (ibnd == 0) { g[0] = U(1, 1); } else { g[0] = U(1, npts); } return; } static void NAG_CALL monit1(double t, Integer npts, Integer npde, const double x[], const double u[], double fmon[], Nag_Comm *comm) { double h1, h2, h3; Integer i; for (i = 2; i <= npts-1; ++i) { h1 = x[i - 1] - x[i - 2]; h2 = x[i] - x[i - 1]; h3 = 0.5* (x[i] - x[i - 2]); /* Second derivatives */ fmon[i-1] = fabs(((U(1, i+1) - U(1, i))/h2 - (U(1, i) - U(1, i-1))/h1)/h3); } fmon[0] = fmon[1]; fmon[npts-1] = fmon[npts-2]; 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] = uleft[0]; return; } int ex2(void) { const Integer npde = 1, npts = 61, ncode = 0, nxi = 0, nxfix = 0, itype = 1; const Integer neqn = npde*npts+ ncode, intpts = 7, lisave = 25+nxfix+neqn; const Integer nwkres = npde*(3*npts+3*npde+32)+7*npts+3, lenode = 11*neqn+50; const Integer mlu = 3*npde-1, lrsave = (3*mlu+1)*neqn+nwkres+lenode; static double xout[7] = { 0., .3, .4, .5, .6, .7, 1. }; double con, dxmesh, tout, trmesh, ts, xratio; Integer exit_status, i, ind, ipminf, it, itask, itol, itrace, m, nrmesh; Nag_Boolean remesh; double *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0, *ue = 0; double *uout = 0, *x = 0, *xfix = 0, *xi = 0; Integer *isave = 0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; INIT_FAIL(fail); exit_status = 0; /* 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(npts, double)) || !(ue = NAG_ALLOC(npde*intpts, double)) || !(uout = NAG_ALLOC(npde*intpts*itype, double)) || !(x = NAG_ALLOC(npts, double)) || !(xfix = NAG_ALLOC(1, double)) || !(xi = NAG_ALLOC(1, double)) || !(isave = NAG_ALLOC(lisave, Integer))) { printf("Allocation failure\n"); exit_status = 1; goto END; } printf("\n\nExample 2\n\n"); itrace = 0; itol = 1; atol [0] = 5e-4; rtol [0] = 0.05; 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); xfix [0] = 0.0; /* Set remesh parameters */ remesh = Nag_TRUE; nrmesh = 5; dxmesh = 0.0; trmesh = 0.0; con = 1.0/(npts-1.0); xratio = 1.5; ipminf = 0; 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.0025; ts = 0.0; tout = 0.0; for (it = 0; it < 2; ++it) { tout = 0.2* (it+1); /* nag_pde_parab_1d_cd_ode_remesh (d03psc), see above. */ nag_pde_parab_1d_cd_ode_remesh(npde, &ts, tout, pdef2, nmflx2, bndry2, uvin2, u, npts, x, ncode, NULLFN, nxi, xi, neqn, rtol, atol, itol, Nag_OneNorm, Nag_LinAlgBand, algopt, remesh, nxfix, xfix, nrmesh, dxmesh, trmesh, ipminf, xratio, con, monit2, 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_remesh (d03psc).\n%s\n", fail.message); exit_status = 1; goto END; } printf(" t = %6.3f\n", ts); printf(" x Approx u Exact u\n\n"); /* Interpolate at output points */ m = 0; /* nag_pde_interp_1d_fd (d03pzc), see above. */ nag_pde_interp_1d_fd(npde, m, u, npts, x, xout, intpts, itype, uout, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_pde_interp_1d_fd (d03pzc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Check against exact solution */ exact(tout, ue, xout, npde, intpts); for (i = 1; i <= intpts; ++i) { printf(" %9.4f", xout[i-1]); printf(" %9.4f", UOUT(1, i, 1)); printf(" %9.4f\n", UE(1, i)); } } 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 (xfix) NAG_FREE(xfix); if (xi) NAG_FREE(xi); if (isave) NAG_FREE(isave); return exit_status; } static void NAG_CALL uvin2(Integer npde, Integer npts, Integer nxi, const double x[], const double xi[], double u[], Integer ncode, double v[], Nag_Comm *comm) { double t; t = 0.0; exact(t, u, x, npde, npts); return; } static void NAG_CALL pdef2(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) { P(1, 1) = 1.0; c[0] = 0.0; d[0] = 0.0; s[0] = -100.0*u[0]*(u[0]-1.0)*(u[0]-0.5); 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) { /* Solution known to be constant at both boundaries */ double ue[1]; if (ibnd == 0) { exact(t, ue, &x[0], npde, 1); g [ 0 ] = UE(1, 1) - U(1, 1); } else { exact(t, ue, &x[npts-1], npde, 1); g [ 0 ] = UE(1, 1) - U(1, npts); } 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) { flux [ 0 ] = uleft[0]; return; } static void NAG_CALL monit2(double t, Integer npts, Integer npde, const double x [], const double u [], double fmon[], Nag_Comm *comm) { static double xa = 0.0; static Integer icount = 0; double h1, ux, uxmax, xl, xleft, xmax, xr, xright; Integer i; /* Locate shock */ uxmax = 0.0; xmax = 0.0; for (i = 2; i <= npts-1; ++i) { h1 = x [i - 1] - x[i - 2]; ux = fabs((U(1, i) - U(1, i - 1))/h1); if (ux > uxmax) { uxmax = ux; xmax = x [i - 1]; } } /* Assign width (on first call only) */ if (icount == 0) { icount = 1; xleft = xmax - x[0]; xright = x [npts-1] - xmax; if (xleft > xright) { xa = xright; } else { xa = xleft; } } xl = xmax - xa; xr = xmax + xa; /* Assign monitor function */ for (i = 0; i < npts; ++i) { if (x [i] > xl && x[i] < xr) { fmon [ i ] = 1.0 + cos(nag_pi *(x[i] - xmax)/xa); } else { fmon [ i ] = 0.0; } } return; } static void NAG_CALL exact(double t, double *u, const double *x, Integer npde, Integer npts) { /* Exact solution (for comparison and b.c. purposes) */ double del, psi, rm, rn, s; Integer i; s = 0.1; del = 0.01; rm = -1.0/del; rn = s /del + 1.0; for (i = 1; i <= npts; ++i) { psi = x[i-1] - t; if (psi < s) { U(1, i) = 1.0; } else if (psi > del + s) { U(1, i) = 0.0; } else { U(1, i) = rm*psi + rn; } } return; }