1. Matthew Emmett
  2. petsc

Commits

Matthew Emmett  committed e9e7197

TSSDC: Tidy/update comments.

  • Participants
  • Parent commits d8e6e14
  • Branches memmett/tssdc

Comments (0)

Files changed (3)

File src/ts/impls/sdc/poly.c

View file
 /*
   Polynomial generation and manipulation routines.
+
+  These polynomial routines are used by the routines in quadrature.c to help
+  generate collocation nodes and SDC integration matrices at runtime based on
+  the number (nnodes) and type (qtype) of SDC nodes requested by -ts_sdc_nodes
+  and -ts_sdc_type.
+
+  Polynomials are stored as simple arrays of coefficients, with the first
+  coeffecient corresponding to the constant term of the polynomial.
 */
 #define PETSC_DESIRE_COMPLEX
 #include <petsc-private/petscimpl.h>                /*I   "petscts.h"   I*/
 
 
 /*
-  Evaluate polynomial.
+  Evaluate polynomial p of degree n at x.
 */
 PetscReal poly_eval(PetscReal *p, PetscInt n, PetscReal x)
 {

File src/ts/impls/sdc/quadrature.c

View file
 /*
   SDC collocation node and quadrature matrix generation routines.
+
+  These polynomial routines generate collocation nodes and SDC integration
+  matrices at runtime based on the number (nnodes) and type (qtype) of SDC nodes
+  requested by -ts_sdc_nodes and -ts_sdc_type.
 */
 #define PETSC_DESIRE_COMPLEX
 #include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/
 
-// XXX: not sure where these belong...
-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);
+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);
 
       if (((flags[m] & SDC_NODE_PROPER) == 0) || (m == i))
         continue;
 
-      // p_{m+1}(x) = (x - x_j) * p_m(x)
+      /* 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];
       for (j=0; j<nsrc+1; j++) p1[j]   -= p[j] * src[m];
       for (j=0; j<nsrc+1; j++) p[j]     = p1[j];
     }
 
-    // evaluate integrals
+    /* evaluate integrals */
     den = poly_eval(p, nsrc, src[i]);
     poly_int(p, nsrc+1);
 
 
   PetscFunctionReturn(0);
 }
-

File src/ts/impls/sdc/sdc.c

View file
 /*
-  Code for timestepping with Spectral Deferred Corredtion IMEX method
+  Code for timestepping with Spectral Deferred Correction (SDC) method.
 
   Notes:
   The general system is written as
 
   where F represents the stiff part of the physics and G represents the non-stiff part.
 
-  XXX: Have assumed that general system is in fact:
+  Note/question from MWE: I have assumed that general system is in fact:
 
     Udot = f^I(t, U) + f^E(t, U)
 
-  so that F(t,U,Udot) = Udot - f^t(t,U) and G(t,U) = f^E(t, U).
+  so that F(t,U,Udot) = Udot - f^t(t,U) and G(t,U) = f^E(t, U).  Is this reasonable?
 
 */
 #include <petsc-private/tsimpl.h>                /*I   "petscts.h"   I*/
 static TSSDCType TSSDCMethodDefault = IMEX;
 
 typedef struct {
-  TSSDCType    qtype;
-  TSSDCMethod  method;
-  PetscInt     nnodes;
-  PetscInt     niters;
-  PetscInt     *flags;
-  PetscReal    *nodes;
+  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 */
 
 
 /*
- Perform one SDC updating sweep.
+ Perform one SDC updating sweep/iteration.
  */
 #undef __FUNCT__
 #define __FUNCT__ "TSSDCSweep"
   for (k=0; k<sdc->niters; k++) {
     ierr = TSSDCIntegrate(ts);CHKERRQ(ierr);
     ierr = TSSDCSweep(ts);CHKERRQ(ierr);
+    /* @MWE: should compute and check sdc residual here, need advice */
   }
   ierr = TSEvaluateStep(ts,0,ts->vec_sol,NULL);CHKERRQ(ierr);
 
-  ts->ptime     += ts->time_step;
+  ts->ptime += ts->time_step;
   ts->steps++;
 
   PetscFunctionReturn(0);
 
 
   PetscFunctionBegin;
-  ierr = PetscFree(sdc->nodes);CHKERRQ(ierr);
   ierr = PetscFree(sdc->flags);CHKERRQ(ierr);
+  ierr = PetscFree(sdc->nodes);CHKERRQ(ierr);
   ierr = MatDestroy(&sdc->Qmat);CHKERRQ(ierr);
   ierr = MatDestroy(&sdc->Smats[0]);CHKERRQ(ierr);
   ierr = MatDestroy(&sdc->Smats[1]);CHKERRQ(ierr);
   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);
 }
 
+
 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);
 
   ierr = VecDuplicate(ts->vec_sol,&sdc->RHS);CHKERRQ(ierr);
   ierr = VecDuplicate(ts->vec_sol,&sdc->Udot);CHKERRQ(ierr);
 
-
-  /*
-  MatView(sdc->Qmat, PETSC_VIEWER_STDOUT_SELF);
-  for (i=0; i<3; i++) {
-    MatView(sdc->Smats[i], PETSC_VIEWER_STDOUT_SELF);
-  }
-   */
-
   PetscFunctionReturn(0);
 }
 
 #undef __FUNCT__
 #define __FUNCT__ "TSSDCGetMethod"
 /*@C
-  TSSDCGetMethod - Get the SDC substepping method (explicit, implicit imex).
+  TSSDCGetMethod - Get the SDC substepping method (explicit, implicit, or imex).
 
   Logically collective
 
 
 .seealso: TSSDCSetMethod()
 @*/
-PetscErrorCode  TSSDCGetMethod(TS ts,TSSDCMethod *sdcmethod)
+PetscErrorCode TSSDCGetMethod(TS ts,TSSDCMethod *sdcmethod)
 {
-  TS_SDC         *sdc = (TS_SDC*)ts->data;
+  TS_SDC *sdc = (TS_SDC*)ts->data;
 
 
   PetscFunctionBegin;
 #undef __FUNCT__
 #define __FUNCT__ "TSSDCSetMethod"
 /*@C
-  TSSDCSetMethod - Set the method of SDC quadrature rule
+  TSSDCSetMethod - Set the SDC substepping method (explicit, implicit, or imex).
 
   Logically collective
 
   Input Parameter:
 +  ts - timestepping context
--  sdcmethod - method of SDC quadrature rule
+-  sdcmethod - method used for substepping
 
   Level: intermediate
 
 .seealso: TSSDCGetMethod(), TSSDCMethod
 @*/
-PetscErrorCode  TSSDCSetMethod(TS ts,TSSDCMethod sdcmethod)
+PetscErrorCode TSSDCSetMethod(TS ts,TSSDCMethod sdcmethod)
 {
-  TS_SDC         *sdc = (TS_SDC*)ts->data;
+  TS_SDC *sdc = (TS_SDC*)ts->data;
 
 
   PetscFunctionBegin;
 
 .seealso: TSSDCSetType()
 @*/
-PetscErrorCode  TSSDCGetType(TS ts,TSSDCType *sdctype)
+PetscErrorCode TSSDCGetType(TS ts,TSSDCType *sdctype)
 {
-  TS_SDC         *sdc = (TS_SDC*)ts->data;
+  TS_SDC *sdc = (TS_SDC*)ts->data;
 
 
   PetscFunctionBegin;
 
 .seealso: TSSDCGetType(), TSSDCType
 @*/
-PetscErrorCode  TSSDCSetType(TS ts,TSSDCType sdctype)
+PetscErrorCode TSSDCSetType(TS ts,TSSDCType sdctype)
 {
-  TS_SDC         *sdc = (TS_SDC*)ts->data;
+  TS_SDC *sdc = (TS_SDC*)ts->data;
 
 
   PetscFunctionBegin;
   ts->ops->load           = TSLoad_SDC;
   ts->ops->setup          = TSSetUp_SDC;
   ts->ops->step           = TSStep_SDC;
-  ts->ops->interpolate    = PETSC_NULL;            /* XXX: this should be easy since we have a collocation solution */
+  ts->ops->interpolate    = PETSC_NULL;            /* @MWE: finish this */
   ts->ops->evaluatestep   = TSEvaluateStep_SDC;
   ts->ops->setfromoptions = TSSetFromOptions_SDC;
   ts->ops->snesfunction   = SNESTSFormFunction_SDC;