/* nag_pde_parab_1d_keller_ode_remesh (d03prc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL pdedef(Integer, double, double, const double[], const double[], const double[], Integer, const double[], const double[], double[], Integer *, Nag_Comm *); static void NAG_CALL bndary(Integer, double, Integer, Integer, const double[], const double[], Integer, const double[], const double[], double[], Integer *, Nag_Comm *); static void NAG_CALL uvinit(Integer, Integer, Integer, const double[], const double[], double[], Integer, double[], Nag_Comm *); static void NAG_CALL monitf(double, Integer, Integer, const double[], const double[], double[], Nag_Comm *); #ifdef __cplusplus } #endif static void NAG_CALL exact(double, Integer, Integer, double *, double *); #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(int argc, char *argv[]) { FILE *fpout; const Integer npde = 2, npts = 61, ncode = 0, nxi = 0, nxfix = 0, nleft = 1; const Integer itype = 1, intpts = 5, neqn = npde*npts+ncode; const Integer lisave = 25+nxfix; const Integer nwkres = npde*(npts+3*npde+21)+7*npts+nxfix+3; const Integer lenode = 11*neqn+50, lrsave = neqn*neqn+neqn+nwkres+lenode; double con, dxmesh, tout, trmesh, ts, xratio; Integer exit_status = 0, i, ind, ipminf, it, itask, itol, itrace, nrmesh; Nag_Boolean remesh, theta; double *algopt = 0, *atol = 0, *rsave = 0, *rtol = 0, *u = 0, *ue = 0; double *uout = 0, *x = 0, *xfix = 0, *xi = 0, *xout = 0; Integer *isave = 0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; INIT_FAIL(fail); /* Check for command-line IO options */ fpout = nag_example_file_io(argc, argv, "-results", NULL); fprintf(fpout, "nag_pde_parab_1d_keller_ode_remesh (d03prc) Example Program" " Results\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*npts, double)) || !(uout = NAG_ALLOC(npde*intpts*itype, double)) || !(x = NAG_ALLOC(npts, double)) || !(xfix = NAG_ALLOC(1, double)) || !(xi = NAG_ALLOC(1, double)) || !(xout = NAG_ALLOC(intpts, double)) || !(isave = NAG_ALLOC(lisave, Integer))) { fprintf(fpout, "Allocation failure\n"); exit_status = 1; goto END; } itrace = 0; itol = 1; atol[0] = 5.0e-5; rtol[0] = atol[0]; fprintf(fpout, " Accuracy requirement =%12.3e", atol[0]); fprintf(fpout, " Number of points = %3ld\n\n", npts); /* Set remesh parameters */ remesh = Nag_TRUE; nrmesh = 3; dxmesh = 0.0; trmesh = 0.0; con = 5.0/(npts-1.0); xratio = 1.2; ipminf = 0; fprintf(fpout, " Remeshing every %3ld time steps\n\n", nrmesh); /* Initialise mesh */ for (i = 0; i < npts; ++i) x[i] = i/(npts-1.0); xout[0] = 0.0; xout[1] = 0.25; xout[2] = 0.5; xout[3] = 0.75; xout[4] = 1.0; fprintf(fpout, " x "); for (i = 0; i < intpts; ++i) { fprintf(fpout, "%10.4f", xout[i]); fprintf(fpout, (i+1)%5 == 0 || i == 4?"\n":" "); } fprintf(fpout, "\n\n"); xi[0] = 0.0; ind = 0; itask = 1; /* Set theta to TRUE if the Theta integrator is required */ theta = Nag_FALSE; for (i = 0; i < 30; ++i) algopt[i] = 0.0; if (theta) { algopt[0] = 2.0; algopt[5] = 2.0; algopt[6] = 1.0; } /* Loop over output value of t */ ts = 0.0; tout = 0.0; for (it = 0; it < 5; ++it) { tout = 0.05*(it+1); /* nag_pde_parab_1d_keller_ode_remesh (d03prc). * General system of first-order PDEs, coupled DAEs, method * of lines, Keller box discretisation, remeshing, one space * variable */ nag_pde_parab_1d_keller_ode_remesh(npde, &ts, tout, pdedef, bndary, uvinit, u, npts, x, nleft, ncode, d03pek, nxi, xi, neqn, rtol, atol, itol, Nag_TwoNorm, Nag_LinAlgFull, algopt, remesh, nxfix, xfix, nrmesh, dxmesh, trmesh, ipminf, xratio, con, monitf, rsave, lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail); if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from " "nag_pde_parab_1d_keller_ode_remesh (d03prc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Interpolate at output points */ /* nag_pde_interp_1d_fd (d03pzc). PDEs, spatial interpolation with * nag_pde_parab_1d_keller_ode_remesh (d03prc). */ nag_pde_interp_1d_fd(npde, 0, u, npts, x, xout, intpts, itype, uout, &fail); if (fail.code != NE_NOERROR) { fprintf(fpout, "Error from nag_pde_interp_1d_fd (d03pzc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Check against exact solution */ exact(ts, npde, intpts, xout, ue); fprintf(fpout, " t = %6.3f\n", ts); fprintf(fpout, " Approx u1"); for (i = 1; i <= intpts; ++i) { fprintf(fpout, "%10.4f", UOUT(1, i, 1)); fprintf(fpout, i%5 == 0 || i == 5?"\n":""); } fprintf(fpout, " Exact u1"); for (i = 1; i <= 5; ++i) { fprintf(fpout, "%10.4f", UE(1, i)); fprintf(fpout, i%5 == 0 || i == 5?"\n":""); } fprintf(fpout, " Approx u2"); for (i = 1; i <= 5; ++i) { fprintf(fpout, "%10.4f", UOUT(2, i, 1)); fprintf(fpout, i%5 == 0 || i == 5?"\n":""); } fprintf(fpout, " Exact u2"); for (i = 1; i <= 5; ++i) { fprintf(fpout, "%10.4f", UE(2, i)); fprintf(fpout, i%5 == 0 || i == 5?"\n":""); } fprintf(fpout, "\n"); } fprintf(fpout, " Number of integration steps in time = %6ld\n", isave[0]); fprintf(fpout, " Number of function evaluations = %6ld\n", isave[1]); fprintf(fpout, " Number of Jacobian evaluations =%6ld\n", isave[2]); fprintf(fpout, " Number of iterations = %6ld\n\n", isave[4]); END: if (fpout != stdout) fclose(fpout); 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 (xout) NAG_FREE(xout); if (isave) NAG_FREE(isave); return exit_status; } static void NAG_CALL uvinit(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) { U(1, i) = exp(x[i-1]); U(2, i) = x[i-1]*x[i-1] + sin(2.0*nag_pi*(x[i-1]*x[i-1])); } return; } static void NAG_CALL pdedef(Integer npde, double t, double x, const double u[], const double udot[], const double ux[], Integer ncode, const double v[], const double vdot[], double res[], Integer *ires, Nag_Comm *comm) { if (*ires == -1) { res[0] = udot[0]; res[1] = udot[1]; } else { res[0] = udot[0] + ux[0] + ux[1]; res[1] = udot[1] + 4.0*ux[0] + ux[1]; } return; } static void NAG_CALL bndary(Integer npde, double t, Integer ibnd, Integer nobc, const double u[], const double udot[], Integer ncode, const double v[], const double vdot[], double res[], Integer *ires, Nag_Comm *comm) { double pp; pp = 2.0*nag_pi; if (ibnd == 0) { if (*ires == -1) { res[0] = 0.0; } else { res[0] = u[0] - 0.5*(exp(t) + exp(-3.0*t)) - 0.25*(sin(9.0*pp*t* t) - sin(pp*t*t)) - 2.0*t*t; } } else { if (*ires == -1) { res[0] = 0.0; } else { res[0] = u[1] - (exp(1.0-3.0*t) - exp(1.0+t) + 0.5*sin(pp*(1.0-3.0*t)*(1.0-3.0*t)) + 0.5*sin(pp*(1.0+t)*(1.0+t)) + 1.0 + 5.0*t*t - 2.0*t); } } return; } static void NAG_CALL monitf(double t, Integer npts, Integer npde, const double x[], const double u[], double fmon[], Nag_Comm *comm) { double d2x1, d2x2, 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 */ d2x1 = fabs(((U(1, i+1)-U(1, i))/h2-(U(1, i)-U(1, i-1))/h1)/h3); d2x2 = fabs(((U(2, i+1)-U(2, i))/h2-(U(2, i)-U(2, i-1))/h1)/h3); fmon[i-1] = d2x1; if (d2x2 > d2x1) fmon[i- 1] = d2x2; } fmon[0] = fmon[1]; fmon[npts- 1] = fmon[npts-2]; return; } static void NAG_CALL exact(double t, Integer npde, Integer npts, double *x, double *u) { /* Exact solution (for comparison purposes) */ double pp; Integer i; pp = 2.0*nag_pi; for (i = 1; i <= npts; ++i) { U(1, i) = 0.5* (exp(x[i-1]+t) + exp(x[i-1]-3.0*t)) + 0.25*(sin(pp*(x[i-1]-3.0*t)*(x[i-1]-3.0*t)) - sin(pp*(x[i-1]+t)*(x[i-1]+t))) + 2.0* t*t - 2.0*x[i-1]*t; U(2, i) = exp(x[i-1]-3.0*t) - exp(x[i-1]+t) + 0.5*(sin(pp*((x[i-1]-3.0*t)*(x[i-1]-3.0*t))) + sin(pp*((x[i-1]+t)*(x[i-1]+t)))) + x[i-1]*x[i-1] + 5.0*t*t - 2.0*x[i-1]*t; } return; }