Commits

Matthew Emmett committed d4e3c23

TSSDC: Add sdc.h, remove lots of whitespace, rename poly routines.

Comments (0)

Files changed (4)

src/ts/impls/sdc/poly.c

 #define PETSC_DESIRE_COMPLEX
 #include <petsc-private/petscimpl.h>                /*I   "petscts.h"   I*/
 #include <petscmath.h>
+#include <../src/ts/impls/sdc/sdc.h>
 
 
 /*
   Evaluate polynomial p of degree n at x.
 */
-PetscReal poly_eval(PetscReal *p, PetscInt n, PetscReal x)
+PetscReal TSSDCPolyEval(PetscReal *p, PetscInt n, PetscReal x)
 {
   PetscInt  j;
   PetscReal v;
-
   v = p[n];
-  for (j=n-1; j>=0; j--)
-    v = x * v + p[j];
-
+  for (j=n-1; j>=0; j--) v = x*v + p[j];
   return v;
 }
 
-
-PetscComplex poly_eval_complex(PetscReal *p, PetscInt n, PetscComplex x)
+PetscComplex TSSDCPolyEvalComplex(PetscReal *p, PetscInt n, PetscComplex x)
 {
   PetscInt     j;
   PetscComplex v;
-
   v = p[n];
-  for (j=n-1; j>=0; j--)
-    v = x * v + p[j];
-
+  for (j=n-1; j>=0; j--) v = x*v + p[j];
   return v;
 }
 
 /*
   Differentiate polynomial (in place).
 */
-void poly_diff(PetscReal *p, PetscInt n)
+void TSSDCPolyDiff(PetscReal *p, PetscInt n)
 {
   PetscInt j;
-  for (j=1; j<n+1; j++)
-    p[j-1] = j * p[j];
+  for (j=1; j<n+1; j++) p[j-1] = j*p[j];
   p[n] = 0.0;
 }
 
 /*
   Integrate polynomial (in place).
  */
-void poly_int(PetscReal *p, PetscInt n)
+void TSSDCPolyInt(PetscReal *p, PetscInt n)
 {
   PetscInt j;
-  for (j=n-2; j>=0; j--)
-    p[j+1] = p[j] / (j+1);
+  for (j=n-2; j>=0; j--) p[j+1] = p[j]/(j+1);
   p[0] = 0.0;
 }
 
   formula.
 */
 #undef __FUNCT__
-#define __FUNCT__ "poly_legendre"
-PetscErrorCode poly_legendre(PetscReal *p, PetscInt n)
+#define __FUNCT__ "TSSDCPolyLegendre"
+PetscErrorCode TSSDCPolyLegendre(PetscReal *p, PetscInt n)
 {
-  PetscInt  j, m, ierr;
-  PetscReal *p0, *p1, *p2;
-
+  PetscInt  j,m,ierr;
+  PetscReal *p0,*p1,*p2;
 
   PetscFunctionBegin;
-
-  ierr = PetscMalloc((n+1)*sizeof(PetscReal), &p0);CHKERRQ(ierr);
-  ierr = PetscMalloc((n+1)*sizeof(PetscReal), &p1);CHKERRQ(ierr);
-  ierr = PetscMalloc((n+1)*sizeof(PetscReal), &p2);CHKERRQ(ierr);
+  ierr = PetscMalloc3(n+1,PetscReal,&p0,n+1,PetscReal,&p1,n+1,PetscReal,&p2);CHKERRQ(ierr);
 
   if (n == 0) {
     p[0] = 1.0;
-    return 0;
+    PetscFunctionReturn(0);
   }
-
   if (n == 1) {
     p[0] = 0.0;
     p[1] = 1.0;
-    return 0;
+    PetscFunctionReturn(0);
   }
 
   /* initialize polynomials, set P_0 = 1, P_1 = x */
-  for (j=0; j<n+1; j++)
-    p0[j] = p1[j] = p2[j] = 0.0;
-
+  for (j=0; j<n+1; j++) p0[j] = p1[j] = p2[j] = 0.0;
   p0[0] = 1.0;
   p1[1] = 1.0;
 
   /* (n + 1) P_{n+1} = (2n + 1) x P_{n} - n P_{n-1} */
   for (m=1; m<n; m++) {
-    for (j=1; j<n+1; j++)
-      p2[j] = ( (2*m + 1) * p1[j-1] - m * p0[j] ) / (m + 1);
+    for (j=1; j<n+1; j++) p2[j] = ( (2*m + 1) * p1[j-1] - m * p0[j] ) / (m + 1);
     p2[0] = - m * p0[0] / (m + 1);
-
     for (j=0; j<n+1; j++) {
       p0[j] = p1[j];
       p1[j] = p2[j];
     }
   }
+  for (j=0; j<n+1; j++) p[j] = p2[j];
 
-  for (j=0; j<n+1; j++)
-    p[j] = p2[j];
-
-  ierr = PetscFree(p0);CHKERRQ(ierr);
-  ierr = PetscFree(p1);CHKERRQ(ierr);
-  ierr = PetscFree(p2);CHKERRQ(ierr);
-
+  ierr = PetscFree3(p0,p1,p2);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   The roots are assumed to be real.
 */
 #undef __FUNCT__
-#define __FUNCT__ "poly_roots"
-PetscErrorCode poly_roots(PetscReal *roots, PetscReal *p0, PetscInt n)
+#define __FUNCT__ "TSSDCPolyRoots"
+PetscErrorCode TSSDCPolyRoots(PetscReal *roots, PetscReal *p0, PetscInt n)
 {
-  PetscComplex num, den, *z0, *z1;
-  PetscReal    *p, acc;
-  PetscInt     i, j, k, ierr;
-
+  PetscComplex num,den,*z0,*z1;
+  PetscReal    *p,acc;
+  PetscInt     i,j,k,ierr;
 
   PetscFunctionBegin;
+  ierr = PetscMalloc3(n+1,PetscReal,&p,n,PetscComplex,&z0,n,PetscComplex,&z1);CHKERRQ(ierr);
 
-  ierr = PetscMalloc((n+1)*sizeof(PetscReal), &p);CHKERRQ(ierr);
-  ierr = PetscMalloc((n)*sizeof(PetscComplex), &z0);CHKERRQ(ierr);
-  ierr = PetscMalloc((n)*sizeof(PetscComplex), &z1);CHKERRQ(ierr);
-
-  /* normalize p0 */
-  for (j=0; j<n+1; j++)
-    p[j] = p0[j] / p0[n];
-
-  /* initial guess */
+  /* normalize p0, initial guess */
+  for (j=0; j<n+1; j++) p[j] = p0[j] / p0[n];
   for (j=0; j<n; j++) {
     z0[j] = PetscPowComplex(0.4 + 0.9 * I, j);
     z1[j] = z0[j];
   /* durand-kerner-weierstrass iterations */
   for (k=0; k<100; k++) {
     for (i=0; i<n; i++) {
-
-      num = poly_eval_complex(p, n, z0[i]);
-
+      num = TSSDCPolyEvalComplex(p, n, z0[i]);
       den = 1.0;
       for (j=0; j<n; j++) {
         if (j == i) continue;
-        den = den * (z0[i] - z0[j]);
+        den *= z0[i] - z0[j];
       }
-
-      z0[i] = z0[i] - num / den;
+      z0[i] -= num / den;
     }
 
     /* converged? */
     acc = 0.0;
-    for (j=0; j<n; j++)
-      acc += PetscAbsComplex(z0[j] - z1[j]);
-
+    for (j=0; j<n; j++) acc += PetscAbsComplex(z0[j] - z1[j]);
     if (acc < PETSC_MACHINE_EPSILON) break;
 
-    for (j=0; j<n; j++)
-      z1[j] = z0[j];
+    for (j=0; j<n; j++) z1[j] = z0[j];
   }
 
-  for (j=0; j<n; j++)
+  for (j=0; j<n; j++) {
     roots[j] = PetscAbsComplex(z0[j]) < PETSC_MACHINE_EPSILON ? 0.0 : PetscRealPart(z0[j]);
-
+  }
   ierr = PetscSortReal(n, roots);CHKERRQ(ierr);
-
-  ierr = PetscFree(p);CHKERRQ(ierr);
-  ierr = PetscFree(z0);CHKERRQ(ierr);
-  ierr = PetscFree(z1);CHKERRQ(ierr);
-
+  ierr = PetscFree3(p,z0,z1);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }

src/ts/impls/sdc/quadrature.c

 */
 #define PETSC_DESIRE_COMPLEX
 #include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/
+#include <../src/ts/impls/sdc/sdc.h>
 
-PetscReal      poly_eval(PetscReal *p, PetscInt n, PetscReal x);
-PetscComplex   poly_eval_complex(PetscReal *p, PetscInt n, PetscComplex x);
-void           poly_diff(PetscReal *p, PetscInt n);
-void           poly_int(PetscReal *p, PetscInt n);
-PetscErrorCode poly_legendre(PetscReal *p, PetscInt n);
-PetscErrorCode poly_roots(PetscReal *roots, PetscReal *p0, PetscInt n);
-
-#define SDC_NODE_PROPER 1
 
 
 #undef __FUNCT__
   PetscInt   j, ierr;
 
   PetscFunctionBegin;
-
   if (qtype == GAUSS_LEGENDRE) {
-
     d = nnodes - 2;
-    ierr = PetscMalloc(sizeof(PetscReal)*(d+1), &p);CHKERRQ(ierr);
-    ierr = PetscMalloc(sizeof(PetscReal)*(d),   &r);CHKERRQ(ierr);
-
-    ierr = poly_legendre(p, d);CHKERRQ(ierr);
-    ierr = poly_roots(r, p, d);CHKERRQ(ierr);
-
+    ierr = PetscMalloc2(d+1,PetscReal,&p,d,PetscReal,&r);CHKERRQ(ierr);
+    ierr = TSSDCPolyLegendre(p, d);CHKERRQ(ierr);
+    ierr = TSSDCPolyRoots(r, p, d);CHKERRQ(ierr);
     nodes[0] = 0.0;
-    for (j=0; j<nnodes-1; j++)
-      nodes[j+1] = 0.5 * (1.0 + r[j]);
+    for (j=0; j<nnodes-1; j++) nodes[j+1] = 0.5 * (1.0 + r[j]);
     nodes[nnodes-1] = 1.0;
-
-    for (j=1; j<nnodes-1; j++)
-      flags[j] |= SDC_NODE_PROPER;
-
-    ierr = PetscFree(r);CHKERRQ(ierr);
-    ierr = PetscFree(p);CHKERRQ(ierr);
-
+    for (j=1; j<nnodes-1; j++) flags[j] |= SDC_NODE_PROPER;
+    ierr = PetscFree2(r,p);CHKERRQ(ierr);
   } else if (qtype == GAUSS_LOBATTO) {
-
     d = nnodes - 1;
-    ierr = PetscMalloc(sizeof(PetscReal)*(d+1), &p);CHKERRQ(ierr);
-    ierr = PetscMalloc(sizeof(PetscReal)*(d-1), &r);CHKERRQ(ierr);
-
-    ierr = poly_legendre(p, d);CHKERRQ(ierr);
-    poly_diff(p, d);
-    ierr = poly_roots(r, p, d-1);CHKERRQ(ierr);
-
+    ierr = PetscMalloc2(d+1,PetscReal,&p,d-1,PetscReal,&r);CHKERRQ(ierr);
+    ierr = TSSDCPolyLegendre(p, d);CHKERRQ(ierr);
+    TSSDCPolyDiff(p, d);
+    ierr = TSSDCPolyRoots(r, p, d-1);CHKERRQ(ierr);
     nodes[0] = 0.0;
-    for (j=0; j<nnodes-2; j++) {
-      nodes[j+1] = 0.5 * (1.0 + r[j]);
-    }
+    for (j=0; j<nnodes-2; j++) nodes[j+1] = 0.5 * (1.0 + r[j]);
     nodes[nnodes-1] = 1.0;
-
-    for (j=0; j<nnodes; j++)
-      flags[j] |= SDC_NODE_PROPER;
-
-    ierr = PetscFree(r);CHKERRQ(ierr);
-    ierr = PetscFree(p);CHKERRQ(ierr);
-
+    for (j=0; j<nnodes; j++) flags[j] |= SDC_NODE_PROPER;
+    ierr = PetscFree2(r,p);CHKERRQ(ierr);
   } else if (qtype == CLENSHAW_CURTIS) {
-
-    for (j=0; j<nnodes; j++)
-      nodes[j] = 0.5 * (1.0 - PetscCosReal(j*PETSC_PI / (nnodes-1)));
-
-    for (j=0; j<nnodes; j++)
-      flags[j] |= SDC_NODE_PROPER;
-
+    for (j=0; j<nnodes; j++) nodes[j] = 0.5 * (1.0 - PetscCosReal(j*PETSC_PI / (nnodes-1)));
+    for (j=0; j<nnodes; j++) flags[j] |= SDC_NODE_PROPER;
   } else {
-
     SETERRQ(0, PETSC_ERR_ARG_WRONG, "Invalid SDC quadrature rule type (qtype)");
-
   }
-
   PetscFunctionReturn(0);
 }
 
 #define __FUNCT__ "TSSDCMats"
 PetscErrorCode TSSDCMats(Mat Qmat, Mat Smat, PetscReal *dst, PetscReal *src, PetscInt *flags, PetscInt ndst, PetscInt nsrc)
 {
-  PetscInt   i, j, m, ierr;
-  PetscReal  den, q, s;
-  PetscReal  *p, *p1;
-
+  PetscInt   i,j,m,ierr;
+  PetscReal  den,q,s;
+  PetscReal  *p,*p1;
 
   PetscFunctionBegin;
-
   ierr = MatZeroEntries(Qmat);CHKERRQ(ierr);
   ierr = MatZeroEntries(Smat);CHKERRQ(ierr);
-
-  ierr = PetscMalloc(sizeof(PetscReal)*(nsrc+1), &p);CHKERRQ(ierr);
-  ierr = PetscMalloc(sizeof(PetscReal)*(nsrc+1), &p1);CHKERRQ(ierr);
+  ierr = PetscMalloc2(nsrc+1,PetscReal,&p,nsrc+1,PetscReal,&p1);CHKERRQ(ierr);
 
   for (i=0; i<nsrc; i++) {
-    if ((flags[i] & SDC_NODE_PROPER) == 0)
-      continue;
+    if ((flags[i] & SDC_NODE_PROPER) == 0) continue;
 
     /* construct interpolating polynomial coefficients */
     p[0] = 1.0;
-    for (j=1; j<nsrc+1; j++)
-      p[j] = 0.0;
-
+    for (j=1; j<nsrc+1; j++) p[j] = 0.0;
     for (m=0; m<nsrc; m++) {
-      if (((flags[m] & SDC_NODE_PROPER) == 0) || (m == i))
-        continue;
-
+      if (((flags[m] & SDC_NODE_PROPER) == 0) || (m == i)) continue;
       /* p_{m+1}(x) = (x - x_j) * p_m(x) */
       p1[0] = 0.0;
       for (j=0; j<nsrc; j++)   p1[j+1]  = p[j];
     }
 
     /* evaluate integrals */
-    den = poly_eval(p, nsrc, src[i]);
-    poly_int(p, nsrc+1);
+    den = TSSDCPolyEval(p, nsrc, src[i]);
+    TSSDCPolyInt(p, nsrc+1);
 
     for (j=1; j<ndst; j++) {
-      q = poly_eval(p, nsrc, dst[j]) - poly_eval(p, nsrc, 0.0);
-      s = poly_eval(p, nsrc, dst[j]) - poly_eval(p, nsrc, dst[j-1]);
-
-      MatSetValue(Qmat, j-1, i, q/den, INSERT_VALUES);
-      MatSetValue(Smat, j-1, i, s/den, INSERT_VALUES);
+      q = TSSDCPolyEval(p,nsrc,dst[j]) - TSSDCPolyEval(p,nsrc,0.0);
+      s = TSSDCPolyEval(p,nsrc,dst[j]) - TSSDCPolyEval(p,nsrc,dst[j-1]);
+      ierr = MatSetValue(Qmat,j-1,i,q/den,INSERT_VALUES);CHKERRQ(ierr);
+      ierr = MatSetValue(Smat,j-1,i,s/den,INSERT_VALUES);CHKERRQ(ierr);
     }
   }
-
+  ierr = PetscFree2(p,p1);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }

src/ts/impls/sdc/sdc.c

   so that F(t,U,Udot) = Udot - f^t(t,U) and G(t,U) = f^E(t, U).  Is this reasonable?
 
 */
+#define PETSC_DESIRE_COMPLEX
 #include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/
+#include <../src/ts/impls/sdc/sdc.h>
 
-const char *const TSSDCTypes[] = {"","GAUSS-LEGENDRE","GAUSS-LOBATTO","CLENSHAW-CURTIS"};
+const char *const TSSDCTypes[]   = {"","GAUSS-LEGENDRE","GAUSS-LOBATTO","CLENSHAW-CURTIS"};
 const char *const TSSDCMethods[] = {"","EXPLICIT","IMPLICIT","IMEX"};
 
-static TSSDCType TSSDCTypeDefault = GAUSS_LOBATTO;
+static TSSDCType TSSDCTypeDefault   = GAUSS_LOBATTO;
 static TSSDCType TSSDCMethodDefault = IMEX;
 
-typedef struct {
-  TSSDCType    qtype;           /* type of collocation nodes */
-  TSSDCMethod  method;          /* type of substepping method to use (fE, bE, f/bE) */
-  PetscInt     nnodes;          /* number of sdc/collcation nodes */
-  PetscInt     niters;          /* number of sdc iterations */
-  PetscInt     *flags;          /* node flags (not really used yet) */
-  PetscReal    *nodes;          /* collocation nodes */
-  Mat          Qmat;            /* 0-to-node integration matrix */
-  Mat          Smats[3];        /* node-to-node integration matrices */
-  Vec          *U;              /* solutions at collocation nodes */
-  Vec          *UdotI;          /* implicit function evaluations */
-  Vec          *UdotRHS;        /* explicit function evaluations */
-  Vec          *S;              /* node-to-node corrections */
-  Vec          RHS;             /* sdc rhs */
-  Vec          Udot;            /* work space */
-  PetscReal    t0;
-  PetscReal    dt;
-} TS_SDC;
-
-
-#undef __FUNCT__
-#define __FUNCT__ "TSEvaluateStep_SDC"
-static PetscErrorCode TSEvaluateStep_SDC(TS ts,PetscInt order,Vec U,PetscBool *done)
-{
-  TS_SDC          *sdc = (TS_SDC*)ts->data;
-  PetscErrorCode  ierr;
-
-
-  PetscFunctionBegin;
-  ierr = VecCopy(sdc->U[sdc->nnodes-1], U);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-
 /*
- Compute node-to-node (Smat) integrals of the collocation solution.
+ Compute node-to-node integrals (using Smat) of the collocation solution.
  */
 #undef __FUNCT__
 #define __FUNCT__ "TSSDCIntegrate"
 static PetscErrorCode TSSDCIntegrate(TS ts)
 {
-  TS_SDC          *sdc = (TS_SDC*)ts->data;
-  PetscInt        n, m;
-  PetscErrorCode  ierr;
-  PetscReal       dt;
-  const PetscScalar *se, *si;
-
+  TS_SDC            *sdc = (TS_SDC*)ts->data;
+  PetscInt          n,m;
+  PetscErrorCode    ierr;
+  PetscReal         dt;
+  const PetscScalar *se,*si;
 
   PetscFunctionBegin;
   dt = ts->time_step;
     ierr = VecZeroEntries(sdc->S[n]);CHKERRQ(ierr);
     ierr = MatGetRow(sdc->Smats[1], n, NULL, NULL, &se);CHKERRQ(ierr);
     ierr = MatGetRow(sdc->Smats[2], n, NULL, NULL, &si);CHKERRQ(ierr);
-
     if (sdc->method == EXPLICIT) {
       for (m=0; m<sdc->nnodes; m++) {
         ierr = VecAXPBY(sdc->S[n],dt*se[m],1.0,sdc->UdotRHS[m]);CHKERRQ(ierr);
         ierr = VecAXPBYPCZ(sdc->S[n],dt*se[m],dt*si[m],1.0,sdc->UdotRHS[m],sdc->UdotI[m]);CHKERRQ(ierr);
       }
     }
-
     ierr = MatRestoreRow(sdc->Smats[2], n, NULL, NULL, &si);CHKERRQ(ierr);
     ierr = MatRestoreRow(sdc->Smats[1], n, NULL, NULL, &se);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 
-
 /*
  Perform one SDC updating sweep/iteration.
  */
 static PetscErrorCode TSSDCSweep(TS ts)
 {
   TS_SDC          *sdc = (TS_SDC*)ts->data;
-  PetscInt        m, its, lits;
+  PetscInt        m,its,lits;
   PetscErrorCode  ierr;
 
-
   PetscFunctionBegin;
   sdc->t0 = ts->ptime;
   for (m=0; m<sdc->nnodes-1; m++) {
     sdc->dt = ts->time_step * (sdc->nodes[m+1] - sdc->nodes[m]);
-
-
     if (sdc->method == EXPLICIT) {
-
       ierr = VecWAXPY(sdc->U[m+1],1.0,sdc->S[m],sdc->U[m]);
       ierr = TSComputeRHSFunction(ts,sdc->t0+sdc->dt,sdc->U[m+1],sdc->UdotRHS[m+1]);CHKERRQ(ierr);
-
     } else if (sdc->method == IMPLICIT) {
-
       ierr = VecWAXPY(sdc->RHS,1.0,sdc->S[m],sdc->U[m]);
-
       ierr = SNESSolve(ts->snes,NULL,sdc->U[m+1]);CHKERRQ(ierr);
       ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr);
       ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr);
       ts->snes_its += its; ts->ksp_its += lits;
-
       ierr = TSComputeIFunction(ts,sdc->t0+sdc->dt,sdc->U[m+1],sdc->Udot,sdc->UdotI[m+1],PETSC_FALSE);CHKERRQ(ierr);
-
     } else {
-
       ierr = VecWAXPY(sdc->RHS,sdc->dt,sdc->UdotRHS[m],sdc->U[m]);
       ierr = VecAXPY(sdc->RHS,1.0,sdc->S[m]);
-
       ierr = SNESSolve(ts->snes,NULL,sdc->U[m+1]);CHKERRQ(ierr);
       ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr);
       ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr);
       ts->snes_its += its; ts->ksp_its += lits;
-
       ierr = TSComputeRHSFunction(ts,sdc->t0+sdc->dt,sdc->U[m+1],sdc->UdotRHS[m+1]);CHKERRQ(ierr);
       ierr = TSComputeIFunction(ts,sdc->t0+sdc->dt,sdc->U[m+1],sdc->Udot,sdc->UdotI[m+1],PETSC_FALSE);CHKERRQ(ierr);
-
     }
-
     sdc->t0 += sdc->dt;
   }
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSStep_SDC"
 static PetscErrorCode TSStep_SDC(TS ts)
 {
   TS_SDC          *sdc = (TS_SDC*)ts->data;
-  PetscInt        m, k;
+  PetscInt        m,k;
   PetscErrorCode  ierr;
 
-
   PetscFunctionBegin;
   ierr = VecZeroEntries(sdc->Udot);CHKERRQ(ierr);
   ierr = VecCopy(ts->vec_sol,sdc->U[0]);CHKERRQ(ierr);
       ierr = VecCopy(sdc->UdotI[0],sdc->UdotI[m]);CHKERRQ(ierr);
     }
   }
-
   /* iterate */
   for (k=0; k<sdc->niters; k++) {
     ierr = TSSDCIntegrate(ts);CHKERRQ(ierr);
 
   ts->ptime += ts->time_step;
   ts->steps++;
-
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "TSEvaluateStep_SDC"
+static PetscErrorCode TSEvaluateStep_SDC(TS ts,PetscInt order,Vec U,PetscBool *done)
+{
+  TS_SDC          *sdc = (TS_SDC*)ts->data;
+  PetscErrorCode  ierr;
+  PetscFunctionBegin;
+  ierr = VecCopy(sdc->U[sdc->nnodes-1], U);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
 
 #undef __FUNCT__
 #define __FUNCT__ "TSReset_SDC"
   TS_SDC         *sdc = (TS_SDC*)ts->data;
   PetscErrorCode ierr;
 
-
   PetscFunctionBegin;
   ierr = PetscFree(sdc->flags);CHKERRQ(ierr);
   ierr = PetscFree(sdc->nodes);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSDestroy_SDC"
 static PetscErrorCode TSDestroy_SDC(TS ts)
 {
   PetscErrorCode ierr;
-
-
   PetscFunctionBegin;
   ierr = TSReset_SDC(ts);CHKERRQ(ierr);
   ierr = PetscFree(ts->data);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
-
 /*
   This defines the nonlinear equation that is to be solved with SNES
   F[t0+dt, U, (U-RHS)/dt] = 0
   TS_SDC         *sdc = (TS_SDC*)ts->data;
   PetscErrorCode ierr;
   PetscReal      invdt = 1./(sdc->dt);
-
-
   PetscFunctionBegin;
   ierr = VecAXPBYPCZ(sdc->Udot,-invdt,invdt,0,sdc->RHS,U);CHKERRQ(ierr);
   ierr = TSComputeIFunction(ts,sdc->t0+sdc->dt,U,sdc->Udot,F,PETSC_FALSE);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "SNESTSFormJacobian_SDC"
 static PetscErrorCode SNESTSFormJacobian_SDC(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *str,TS ts)
   TS_SDC         *sdc = (TS_SDC*)ts->data;
   PetscErrorCode ierr;
   PetscReal      invdt = 1./(sdc->dt);
-
-
   PetscFunctionBegin;
   ierr = TSComputeIJacobian(ts,sdc->t0,U,sdc->Udot,invdt,A,B,str,PETSC_FALSE);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
-
-PetscErrorCode TSSDCNodes(PetscReal *nodes, PetscInt *flags, TSSDCType qtype, PetscInt nnodes);
-PetscErrorCode TSSDCMats(Mat Qmat, Mat Smat, PetscReal *dst, PetscReal *src, PetscInt *flags, PetscInt ndst, PetscInt nsrc);
-
 #undef __FUNCT__
 #define __FUNCT__ "TSSetUp_SDC"
 static PetscErrorCode TSSetUp_SDC(TS ts)
   PetscInt       i;
   PetscErrorCode ierr;
 
-
   PetscFunctionBegin;
 
   /* collocation nodes */
   MatAssemblyBegin(sdc->Smats[2],MAT_FINAL_ASSEMBLY);
   MatAssemblyEnd(sdc->Smats[2],MAT_FINAL_ASSEMBLY);
 
-  /* allocate U, UdotI, UdotRHS */
   ierr = VecDuplicateVecs(ts->vec_sol,sdc->nnodes,&sdc->U);CHKERRQ(ierr);
   ierr = VecDuplicateVecs(ts->vec_sol,sdc->nnodes,&sdc->UdotI);CHKERRQ(ierr);
   ierr = VecDuplicateVecs(ts->vec_sol,sdc->nnodes,&sdc->UdotRHS);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSSetFromOptions_SDC"
 static PetscErrorCode TSSetFromOptions_SDC(TS ts)
   TS_SDC         *sdc       = (TS_SDC*)ts->data;
   PetscErrorCode ierr;
 
-
   PetscFunctionBegin;
   ierr = PetscOptionsHead("SDC ODE solver options");CHKERRQ(ierr);
   {
     PetscInt   choice;
     PetscBool  flg;
-
     ierr = PetscOptionsEList("-ts_sdc_method","Type of SDC method","TSSDCSetMehtod",TSSDCMethods,4,TSSDCMethods[TSSDCMethodDefault],&choice,&flg);CHKERRQ(ierr);
     ierr = TSSDCSetMethod(ts,flg ? (TSSDCMethod)choice : TSSDCMethodDefault);CHKERRQ(ierr);
     ierr = PetscOptionsEList("-ts_sdc_type","Type of SDC quadrature rule","TSSDCSetType",TSSDCTypes,4,TSSDCTypes[TSSDCTypeDefault],&choice,&flg);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSView_SDC"
 static PetscErrorCode TSView_SDC(TS ts,PetscViewer viewer)
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSLoad_SDC"
 static PetscErrorCode TSLoad_SDC(TS ts,PetscViewer viewer)
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSSDCGetMethod"
 /*@C
 PetscErrorCode TSSDCGetMethod(TS ts,TSSDCMethod *sdcmethod)
 {
   TS_SDC *sdc = (TS_SDC*)ts->data;
-
-
   PetscFunctionBegin;
   *sdcmethod = sdc->method;
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSSDCSetMethod"
 /*@C
 PetscErrorCode TSSDCSetMethod(TS ts,TSSDCMethod sdcmethod)
 {
   TS_SDC *sdc = (TS_SDC*)ts->data;
-
-
   PetscFunctionBegin;
   sdc->method = sdcmethod;
   PetscFunctionReturn(0);
 PetscErrorCode TSSDCGetType(TS ts,TSSDCType *sdctype)
 {
   TS_SDC *sdc = (TS_SDC*)ts->data;
-
-
   PetscFunctionBegin;
   *sdctype = sdc->qtype;
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "TSSDCSetType"
 /*@C
 PetscErrorCode TSSDCSetType(TS ts,TSSDCType sdctype)
 {
   TS_SDC *sdc = (TS_SDC*)ts->data;
-
-
   PetscFunctionBegin;
   sdc->qtype = sdctype;
   PetscFunctionReturn(0);
 }
 
-
 /*MC
       TSSDC - ODE solver using Spectral Deferred Correction schemes
 
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-
   ts->ops->reset          = TSReset_SDC;
   ts->ops->destroy        = TSDestroy_SDC;
   ts->ops->view           = TSView_SDC;
   ts->ops->setfromoptions = TSSetFromOptions_SDC;
   ts->ops->snesfunction   = SNESTSFormFunction_SDC;
   ts->ops->snesjacobian   = SNESTSFormJacobian_SDC;
-
   ierr = PetscNewLog(ts,TS_SDC,&th);CHKERRQ(ierr);
   ts->data = (void*)th;
   th->nnodes = 3;
   th->niters = 4;
-
   PetscFunctionReturn(0);
 }

src/ts/impls/sdc/sdc.h

+#ifndef __PETSCSDC_H
+#define __PETSCSDC_H
+
+#define PETSC_DESIRE_COMPLEX
+#include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/
+#include <petsc-private/petscimpl.h>                /*I   "petscts.h"   I*/
+
+#define SDC_NODE_PROPER (1 << 0)
+
+PetscReal TSSDCPolyEval(PetscReal *p, PetscInt n, PetscReal x);
+PetscComplex TSSDCPolyEvalComplex(PetscReal *p, PetscInt n, PetscComplex x);
+void TSSDCPolyDiff(PetscReal *p, PetscInt n);
+void TSSDCPolyInt(PetscReal *p, PetscInt n);
+PetscErrorCode TSSDCPolyLegendre(PetscReal *p, PetscInt n);
+PetscErrorCode TSSDCPolyRoots(PetscReal *roots, PetscReal *p0, PetscInt n);
+
+PetscErrorCode TSSDCNodes(PetscReal *nodes, PetscInt *flags, TSSDCType qtype, PetscInt nnodes);
+PetscErrorCode TSSDCMats(Mat Qmat, Mat Smat, PetscReal *dst, PetscReal *src, PetscInt *flags, PetscInt ndst, PetscInt nsrc);
+
+typedef struct {
+  TSSDCType    qtype;           /* type of collocation nodes */
+  TSSDCMethod  method;          /* type of substepping method to use (fE, bE, f/bE) */
+  PetscInt     nnodes;          /* number of sdc/collcation nodes */
+  PetscInt     niters;          /* number of sdc iterations */
+  PetscInt     *flags;          /* node flags (not really used yet) */
+  PetscReal    *nodes;          /* collocation nodes */
+  Mat          Qmat;            /* 0-to-node integration matrix */
+  Mat          Smats[3];        /* node-to-node integration matrices */
+  Vec          *U;              /* solutions at collocation nodes */
+  Vec          *UdotI;          /* implicit function evaluations */
+  Vec          *UdotRHS;        /* explicit function evaluations */
+  Vec          *S;              /* node-to-node corrections */
+  Vec          RHS;             /* sdc rhs */
+  Vec          Udot;            /* work space */
+  PetscReal    t0;
+  PetscReal    dt;
+} TS_SDC;
+
+#endif