# Commits

committed e75a5fb

TS ex33: start segmental refinement example

This is a simplified version of ts/ex11 that solves a Lapacian in 2D
with a 5-point stencil. It has a synthetic solution and a convergence
test. It uses SNES to construct a differencing Jacoabian without a
solve. This is intened to be the start of an extreme scale solver
driver.

• Participants
• Parent commits 4dba586
• Branches jed/next-seaice, jed/sr-driver4 4

# File src/ts/examples/tutorials/ex33.c

+static char help[] = "Second Order Finite Volume Example of Laplacian div \alpha grad u = f.\n";
+/*F
+
+  Multigrid Poisson solver:
+     Laplace(u) = f
+
+  A numerical test on page 64 of
+    A MULTIGRID TUTORIAL by Briggs, Henson & McCormick'
+  u := (x^2-x^4)*(y^4-y^2)*(z^4-z^2) with alpha(x) = -1
+
+  Course grid mesh is read in from an ExodusII file, usually generated by Cubit.
+F*/
+#include <petscts.h>
+#include <petscdmplex.h>
+#include <petscsf.h>
+#include <petscblaslapack.h>
+#if defined(PETSC_HAVE_EXODUSII)
+#include <exodusII.h>
+#else
+#error This example requires ExodusII support. Reconfigure using --download-exodusii
+#endif
+
+#define DIM 2                   /* Geometric dimension */
+#define ALEN(a) (sizeof(a)/sizeof((a)[0]))
+
+/* Represents continuum physical equations. */
+typedef struct _n_Physics *Physics;
+
+/* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
+ * discretization-independent, but its members depend on the scenario being solved. */
+typedef struct _n_Model *Model;
+
+/* 'User' implements a discretization of a continuous model. */
+typedef struct _n_User *User;
+
+typedef PetscErrorCode (*SolutionFunction)(Model,PetscReal,const PetscReal*,PetscScalar*,void*);
+typedef PetscErrorCode (*BoundaryFunction)(Model,PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*);
+typedef PetscErrorCode (*FunctionalFunction)(Model,PetscReal,const PetscReal*,const PetscScalar*,PetscReal*,void*);
+static PetscErrorCode ModelBoundaryRegister(Model,const char*,BoundaryFunction,void*,PetscInt,const PetscInt*);
+static PetscErrorCode ModelSolutionSetDefault(Model,SolutionFunction,void*);
+static PetscErrorCode ModelSourceSetDefault(Model,SolutionFunction,void*);
+static PetscErrorCode ModelAlphaSetDefault(Model,SolutionFunction,void*);
+static PetscErrorCode ModelFunctionalRegister(Model,const char*,PetscInt*,FunctionalFunction,void*);
+static PetscErrorCode OutputVTK(DM,const char*,PetscViewer*);
+
+struct FieldDescription {
+  const char *name;
+  PetscInt dof;
+};
+
+typedef struct _n_BoundaryLink *BoundaryLink;
+struct _n_BoundaryLink {
+  char             *name;
+  BoundaryFunction func;
+  void             *ctx;
+  PetscInt         numids;
+  PetscInt         *ids;
+  BoundaryLink     next;
+};
+
+typedef struct _n_FunctionalLink *FunctionalLink;
+struct _n_FunctionalLink {
+  char               *name;
+  FunctionalFunction func;
+  void               *ctx;
+  PetscInt           offset;
+  FunctionalLink     next;
+};
+
+struct _n_Physics {
+  PetscInt        dof;          /* number of degrees of freedom per cell */
+  void            *data;
+  PetscInt        nfields;
+  const struct FieldDescription *field_desc;
+};
+
+struct _n_Model {
+  MPI_Comm         comm;        /* Does not do collective communicaton, but some error conditions can be collective */
+  Physics          physics;
+  BoundaryLink     boundary;
+  FunctionalLink   functionalRegistry;
+  PetscInt         maxComputed;
+  PetscInt         numMonitored;
+  FunctionalLink   *functionalMonitored;
+  PetscInt         numCall;
+  FunctionalLink   *functionalCall;
+  SolutionFunction solution;
+  SolutionFunction source;
+  SolutionFunction alpha;
+  void             *solutionctx;
+};
+
+struct _n_User {
+  PetscInt       numGhostCells;
+  PetscInt       cEndInterior; /* First boundary ghost cell */
+  Vec            cellgeom, facegeom;
+  PetscInt       vtkInterval;   /* For monitor */
+  Model          model;
+  PetscReal      errorInf[32], error2[32];
+  PetscInt       N[32],ilev;
+};
+
+typedef struct {
+  PetscScalar normal[DIM];              /* Area-scaled normals */
+  PetscScalar centroid[DIM];            /* Location of centroid (quadrature point) */
+} FaceGeom;
+
+typedef struct {
+  PetscScalar centroid[DIM];
+  PetscScalar volume;
+} CellGeom;
+
+PETSC_STATIC_INLINE PetscScalar DotD(const PetscScalar *x, const PetscScalar *y) {PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < DIM; ++d) sum += x[d]*y[d]; return sum;}
+PETSC_STATIC_INLINE PetscReal Norm2D(const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(DotD(x,x)));}
+PETSC_STATIC_INLINE void axD(const PetscScalar a,PetscScalar *x){PetscInt i; for (i=0; i<DIM; i++) x[i] *= a;}
+PETSC_STATIC_INLINE void waxD(const PetscScalar a,const PetscScalar *x, PetscScalar *w){PetscInt i; for (i=0; i<DIM; i++) w[i] = x[i]*a;}
+PETSC_STATIC_INLINE void NormalizeD(PetscScalar *x) {PetscReal a = 1./Norm2D(x); PetscInt d; for (d = 0; d < DIM; ++d) x[d] *= a;}
+PETSC_STATIC_INLINE void WaxpyD(PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w) {PetscInt d; for (d = 0; d < DIM; ++d) w[d] = a*x[d] + y[d];}
+
+/******************* Laplacian ********************/
+typedef struct {
+  struct {
+    PetscInt Error;
+  } functional;
+} Physics_Lap;
+
+static const struct FieldDescription PhysicsFields_Lap[] = {{"U",1},{NULL,0}};
+
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsBoundary_DiriHomo"
+static PetscErrorCode PhysicsBoundary_DiriHomo(Model mod, PetscReal time, const PetscReal *fCentroid, const PetscReal *normal, const PetscScalar *xInterior, PetscScalar *xGhost, void *ctx)
+{
+  /* Physics      phys = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscFunctionBeginUser;
+  xGhost[0] = -xInterior[0];
+  PetscFunctionReturn(0);
+}
+
+/* U: (x^2-x^4)*(y^4-y^2)*(z^4-z^2) */
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsSolution_Lap_x4_x2"
+static PetscErrorCode PhysicsSolution_Lap_x4_x2(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
+{
+  /* Physics        phys    = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscScalar r;
+  PetscScalar e[DIM];
+  PetscInt d;
+
+  PetscFunctionBeginUser;
+  for (d = 0; d < DIM; ++d) e[d] = x[d]*x[d];
+  for (d = 0; d < DIM; ++d) e[d] *= 1 - e[d];
+  for (r = 1.0, d = 0; d < DIM; ++d) r *= e[d];
+  u[0] = r;
+  PetscFunctionReturn(0);
+}
+
+/* RHS: Lap((x^2-x^4)*(y^4-y^2)*(z^4-z^2)) */
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsSource_Lap_x4_x2"
+static PetscErrorCode PhysicsSource_Lap_x4_x2(Model mod,PetscReal time,const PetscReal *x,PetscScalar *rhs,void *ctx)
+{
+  /* Physics        phys    = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscScalar r;
+  PetscScalar x2[DIM]; /* = x*x; */
+  PetscScalar a[DIM]; /* = x2*(1-x2) */
+  PetscScalar b[DIM]; /* = 2-12*x2; */
+  PetscInt d;
+
+  PetscFunctionBeginUser;
+  for (d = 0; d < DIM; ++d) x2[d] = x[d]*x[d];
+  for (d = 0; d < DIM; ++d) a[d] = x2[d]*(1-x2[d]);
+  for (d = 0; d < DIM; ++d) b[d] = 2-12*x2[d];
+
+  /* cross products of coordinates prevents the use of D_TERM */
+  if (DIM==2)
+    r = b[0]*a[1] + a[0]*b[1];
+  else if (DIM==3)
+    r = b[0]*a[1]*a[2] + a[0]*b[1]*a[2] + a[0]*a[1]*b[2];
+  else  if (DIM==1)
+   r = b[0];
+  else
+    SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution this DIM %d",DIM);
+  rhs[0] = -r; /* keep solution positive */
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsAlpha_n1"
+static PetscErrorCode PhysicsAlpha_n1(Model mod,PetscReal time,const PetscReal *x,PetscScalar *alpha,void *ctx)
+{
+  /* Physics        phys    = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscFunctionBeginUser;
+  alpha[0] = -1.0;
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsFunctional_Error"
+static PetscErrorCode PhysicsFunctional_Error(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *u,PetscReal *f,void *ctx)
+{
+  Physics        phys    = (Physics)ctx;
+  Physics_Lap   *lap = (Physics_Lap*)phys->data;
+  PetscScalar    yexact[1];
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = (*mod->solution)(mod,time,x,yexact,phys);CHKERRQ(ierr);
+  /* f[lap->functional.Error] = PetscAbsScalar(u[0]-yexact[0]); */
+  f[lap->functional.Error] = PetscRealPart(u[0]-yexact[0]); /* lets do real error for plotting */
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsCreate_Lap"
+static PetscErrorCode PhysicsCreate_Lap(Model mod,Physics phys)
+{
+  Physics_Lap *lap;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  phys->field_desc = PhysicsFields_Lap;
+  ierr = PetscNew(Physics_Lap,&phys->data);CHKERRQ(ierr);
+  lap = (Physics_Lap*)phys->data;
+
+  {
+    const PetscInt diriids[] = {100,101,102,103};
+    /* Register "canned" boundary conditions and defaults for where to apply. */
+    ierr = ModelBoundaryRegister(mod,"diri",PhysicsBoundary_DiriHomo,phys,ALEN(diriids),diriids);CHKERRQ(ierr);
+    /* Initial/transient solution with default boundary conditions */
+    ierr = ModelSolutionSetDefault(mod,PhysicsSolution_Lap_x4_x2,phys);CHKERRQ(ierr);
+    ierr = ModelSourceSetDefault(mod,PhysicsSource_Lap_x4_x2,phys);CHKERRQ(ierr);
+    ierr = ModelAlphaSetDefault(mod,PhysicsAlpha_n1,phys);CHKERRQ(ierr);
+    /* Register "canned" functionals */
+    ierr = ModelFunctionalRegister(mod,"error",&lap->functional.Error,PhysicsFunctional_Error,phys);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+/* #undef __FUNCT__ */
+/* #define __FUNCT__ "IsExteriorOrGhostFace" */
+/* static PetscErrorCode IsExteriorOrGhostFace(DM dm,PetscInt face,PetscBool *isghost) */
+/* { */
+/*   PetscErrorCode ierr; */
+/*   PetscInt       ghost,boundary; */
+
+/*   PetscFunctionBeginUser; */
+/*   *isghost = PETSC_FALSE; */
+/*   ierr = DMPlexGetLabelValue(dm, "ghost", face, &ghost);CHKERRQ(ierr); */
+/*   ierr = DMPlexGetLabelValue(dm, "Face Sets", face, &boundary);CHKERRQ(ierr); */
+/*   if (ghost >= 0 || boundary >= 0) *isghost = PETSC_TRUE; */
+/*   PetscFunctionReturn(0); */
+/* } */
+
+#undef __FUNCT__
+#define __FUNCT__ "ConstructGeometry"
+/* Set up face data and cell data */
+PetscErrorCode ConstructGeometry(DM dm, Vec *facegeom, Vec *cellgeom, User user)
+{
+  DM             dmFace, dmCell;
+  DMLabel        ghostLabel;
+  PetscSection   sectionFace, sectionCell;
+  PetscSection   coordSection;
+  Vec            coordinates;
+  PetscScalar    *fgeom, *cgeom;
+  PetscInt       dim, cStart, cEnd, c, fStart, fEnd, f;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  if (dim != DIM) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for dim %D != DIM %D",dim,DIM);
+  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+
+  /* Make cell centroids and volumes */
+  ierr = DMClone(dm, &dmCell);CHKERRQ(ierr);
+  ierr = DMPlexSetCoordinateSection(dmCell, coordSection);CHKERRQ(ierr);
+  ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    ierr = PetscSectionSetDof(sectionCell, c, sizeof(CellGeom)/sizeof(PetscScalar));CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&sectionCell);CHKERRQ(ierr); /* relinquish our reference */
+
+  ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr);
+  ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr);
+  for (c = cStart; c < user->cEndInterior; ++c) {
+    CellGeom *cg;
+
+    ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
+    ierr = PetscMemzero(cg,sizeof(*cg));CHKERRQ(ierr);
+    ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr);
+  }
+  /* Compute face normals, centroid(used?) */
+  ierr = DMClone(dm, &dmFace);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionFace);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr);
+
+  for (f = fStart; f < fEnd; ++f) {
+    ierr = PetscSectionSetDof(sectionFace, f, sizeof(FaceGeom)/sizeof(PetscScalar));CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&sectionFace);CHKERRQ(ierr);
+  ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr);
+  ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr);
+  ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr);
+  for (f = fStart; f < fEnd; ++f) {
+    FaceGeom *fg;
+    PetscInt  ghost;
+
+    ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);
+    if (ghost >= 0) continue;
+    ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr);
+    ierr = DMPlexComputeCellGeometryFVM(dm, f, NULL, fg->centroid, fg->normal);CHKERRQ(ierr);
+    /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
+    {
+      CellGeom       *cL, *cR;
+      const PetscInt *cells;
+      PetscReal      *lcentroid, *rcentroid;
+      PetscScalar     v[3];
+      PetscInt        d;
+
+      ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr);
+      ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr);
+      ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr);
+      lcentroid = cells[0] >= user->cEndInterior ? fg->centroid : cL->centroid;
+      rcentroid = cells[1] >= user->cEndInterior ? fg->centroid : cR->centroid;
+      WaxpyD(-1, lcentroid, rcentroid, v);
+      if (DotD(fg->normal, v) < 0) {
+        for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
+      }
+      if (DotD(fg->normal, v) <= 0) {
+#if DIM == 2
+        SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, fg->normal[0], fg->normal[1], v[0], v[1]);
+#elif DIM == 3
+        SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, fg->normal[0], fg->normal[1], fg->normal[2], v[0], v[1], v[2]);
+#else
+#  error DIM not supported
+#endif
+      }
+    }
+  }
+  /* Compute centroids of ghost cells */
+  for (c = user->cEndInterior; c < cEnd; ++c) {
+    FaceGeom       *fg;
+    const PetscInt *cone,    *support;
+    PetscInt        coneSize, supportSize, s;
+
+    ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr);
+    if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
+    ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr);
+    ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr);
+    if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize);
+    ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr);
+    ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr);
+    for (s = 0; s < 2; ++s) {
+      /* Reflect ghost centroid across plane of face */
+      if (support[s] == c) {
+        const CellGeom *ci;
+        CellGeom       *cg;
+        PetscScalar     c2f[3], a;
+
+        ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr);
+        WaxpyD(-1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
+        a    = DotD(c2f, fg->normal)/DotD(fg->normal, fg->normal);
+        ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr);
+        WaxpyD(2*a, fg->normal, ci->centroid, cg->centroid);
+        cg->volume = ci->volume;
+      }
+    }
+  }
+  ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr);
+  ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr);
+  ierr = DMDestroy(&dmCell);CHKERRQ(ierr);
+  ierr = DMDestroy(&dmFace);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SetUpLocalSpace"
+PetscErrorCode SetUpLocalSpace(DM dm, User user)
+{
+  PetscSection   stateSection;
+  Physics        phys;
+  PetscInt       dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, c, i;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection);CHKERRQ(ierr);
+  phys = user->model->physics;
+  ierr = PetscSectionSetNumFields(stateSection,phys->nfields);CHKERRQ(ierr);
+  for (i=0; i<phys->nfields; i++) {
+    ierr = PetscSectionSetFieldName(stateSection,i,phys->field_desc[i].name);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldComponents(stateSection,i,phys->field_desc[i].dof);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetChart(stateSection, cStart, cEnd);CHKERRQ(ierr);
+
+  for (c = cStart; c < cEnd; ++c) {
+    for (i=0; i<phys->nfields; i++) {
+      ierr = PetscSectionSetFieldDof(stateSection,c,i,phys->field_desc[i].dof);CHKERRQ(ierr);
+    }
+    ierr = PetscSectionSetDof(stateSection, c, dof);CHKERRQ(ierr);
+  }
+  for (c = user->cEndInterior; c < cEnd; ++c) {
+    ierr = PetscSectionSetConstraintDof(stateSection, c, dof);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(stateSection);CHKERRQ(ierr);
+  ierr = PetscMalloc(dof * sizeof(PetscInt), &cind);CHKERRQ(ierr);
+  for (d = 0; d < dof; ++d) cind[d] = d;
+  for (c = user->cEndInterior; c < cEnd; ++c) {
+    ierr = PetscSectionSetConstraintIndices(stateSection, c, cind);CHKERRQ(ierr);
+  }
+  ierr = PetscFree(cind);CHKERRQ(ierr);
+  ierr = PetscSectionGetStorageSize(stateSection, &stateSize);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dm,stateSection);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&stateSection);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SetUpBoundaries"
+PetscErrorCode SetUpBoundaries(MPI_Comm comm, User user)
+{
+  Model          mod = user->model;
+  PetscErrorCode ierr;
+  BoundaryLink   b;
+
+  PetscFunctionBeginUser;
+  ierr = PetscOptionsBegin(comm,NULL,"Boundary condition options","");CHKERRQ(ierr);
+  for (b = mod->boundary; b; b=b->next) {
+    char      optname[512];
+    PetscInt  ids[512],len = 512;
+    PetscBool flg;
+    ierr = PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name);CHKERRQ(ierr);
+    ierr = PetscMemzero(ids,sizeof(ids));CHKERRQ(ierr);
+    ierr = PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg);CHKERRQ(ierr);
+    if (flg) {
+      /* TODO: check all IDs to make sure they exist in the mesh */
+      ierr      = PetscFree(b->ids);CHKERRQ(ierr);
+      b->numids = len;
+      ierr      = PetscMalloc(len*sizeof(PetscInt),&b->ids);CHKERRQ(ierr);
+      ierr      = PetscMemcpy(b->ids,ids,len*sizeof(PetscInt));CHKERRQ(ierr);
+    }
+  }
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelBoundaryRegister"
+/* The ids are just defaults, can be overridden at command line. I expect to set defaults based on names in the future. */
+static PetscErrorCode ModelBoundaryRegister(Model mod,const char *name,BoundaryFunction bcFunc,void *ctx,PetscInt numids,const PetscInt *ids)
+{
+  PetscErrorCode ierr;
+  BoundaryLink   link;
+
+  PetscFunctionBeginUser;
+  ierr          = PetscNew(struct _n_BoundaryLink,&link);CHKERRQ(ierr);
+  ierr          = PetscStrallocpy(name,&link->name);CHKERRQ(ierr);
+  link->numids  = numids;
+  ierr          = PetscMalloc(numids*sizeof(PetscInt),&link->ids);CHKERRQ(ierr);
+  ierr          = PetscMemcpy(link->ids,ids,numids*sizeof(PetscInt));CHKERRQ(ierr);
+  link->func    = bcFunc;
+  link->ctx     = ctx;
+  link->next    = mod->boundary;
+  mod->boundary = link;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "BoundaryLinkDestroy"
+static PetscErrorCode BoundaryLinkDestroy(BoundaryLink *link)
+{
+  PetscErrorCode ierr;
+  BoundaryLink   l,next;
+
+  PetscFunctionBeginUser;
+  if (!link) PetscFunctionReturn(0);
+  l     = *link;
+  *link = NULL;
+  for (; l; l=next) {
+    next = l->next;
+    ierr = PetscFree(l->ids);CHKERRQ(ierr);
+    ierr = PetscFree(l->name);CHKERRQ(ierr);
+    ierr = PetscFree(l);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelBoundaryFind"
+static PetscErrorCode ModelBoundaryFind(Model mod,PetscInt id,BoundaryFunction *bcFunc,void **ctx)
+{
+  BoundaryLink link;
+  PetscInt     i;
+
+  PetscFunctionBeginUser;
+  *bcFunc = NULL;
+  for (link=mod->boundary; link; link=link->next) {
+    for (i=0; i<link->numids; i++) {
+      if (link->ids[i] == id) {
+        *bcFunc = link->func;
+        *ctx    = link->ctx;
+        PetscFunctionReturn(0);
+      }
+    }
+  }
+  SETERRQ1(mod->comm,PETSC_ERR_USER,"Boundary ID %D not associated with any registered boundary condition",id);
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "ModelSolutionSetDefault"
+/* Behavior will be different for multi-physics or when using non-default boundary conditions */
+static PetscErrorCode ModelSolutionSetDefault(Model mod,SolutionFunction func,void *ctx)
+{
+  PetscFunctionBeginUser;
+  mod->solution    = func;
+  mod->solutionctx = ctx;
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "ModelSourceSetDefault"
+/* Behavior will be different for multi-physics or when using non-default boundary conditions */
+static PetscErrorCode ModelSourceSetDefault(Model mod,SolutionFunction func,void *ctx)
+{
+  PetscFunctionBeginUser;
+  mod->source    = func;
+  /* mod->solutionctx = ctx; */
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "ModelAlphaSetDefault"
+/* Behavior will be different for multi-physics or when using non-default boundary conditions */
+static PetscErrorCode ModelAlphaSetDefault(Model mod,SolutionFunction func,void *ctx)
+{
+  PetscFunctionBeginUser;
+  mod->alpha    = func;
+  /* mod->solutionctx = ctx; */
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelFunctionalRegister"
+static PetscErrorCode ModelFunctionalRegister(Model mod,const char *name,PetscInt *offset,FunctionalFunction func,void *ctx)
+{
+  PetscErrorCode ierr;
+  FunctionalLink link,*ptr;
+  PetscInt       lastoffset = -1;
+
+  PetscFunctionBeginUser;
+  for (ptr=&mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
+  ierr         = PetscNew(struct _n_FunctionalLink,&link);CHKERRQ(ierr);
+  ierr         = PetscStrallocpy(name,&link->name);CHKERRQ(ierr);
+  link->offset = lastoffset + 1;
+  link->func   = func;
+  link->ctx    = ctx;
+  link->next   = NULL;
+  *ptr         = link;
+  *offset      = link->offset;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelFunctionalSetFromOptions"
+static PetscErrorCode ModelFunctionalSetFromOptions(Model mod)
+{
+  PetscErrorCode ierr;
+  PetscInt       i,j;
+  FunctionalLink link;
+  char           *names[256];
+
+  PetscFunctionBeginUser;
+  mod->numMonitored = ALEN(names);
+  ierr = PetscOptionsStringArray("-monitor","list of functionals to monitor","",names,&mod->numMonitored,NULL);CHKERRQ(ierr);
+  /* Create list of functionals that will be computed somehow */
+  ierr = PetscMalloc(mod->numMonitored*sizeof(FunctionalLink),&mod->functionalMonitored);CHKERRQ(ierr);
+  /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
+  ierr = PetscMalloc(mod->numMonitored*sizeof(FunctionalLink),&mod->functionalCall);CHKERRQ(ierr);
+  mod->numCall = 0;
+
+  for (i=0; i<mod->numMonitored; i++) {
+    for (link=mod->functionalRegistry; link; link=link->next) {
+      PetscBool match;
+      ierr = PetscStrcasecmp(names[i],link->name,&match);CHKERRQ(ierr);
+      if (match) break;
+    }
+    if (!link) SETERRQ1(mod->comm,PETSC_ERR_USER,"No known functional '%s'",names[i]);
+    mod->functionalMonitored[i] = link;
+    for (j=0; j<i; j++) {
+      if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
+    }
+    mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
+next_name:
+    ierr = PetscFree(names[i]);CHKERRQ(ierr);
+  }
+
+  /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
+  mod->maxComputed = -1;
+  for (link=mod->functionalRegistry; link; link=link->next) {
+    for (i=0; i<mod->numCall; i++) {
+      FunctionalLink call = mod->functionalCall[i];
+      if (link->func == call->func && link->ctx == call->ctx) {
+        mod->maxComputed = PetscMax(mod->maxComputed,link->offset);
+      }
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FunctionalLinkDestroy"
+static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
+{
+  PetscErrorCode ierr;
+  FunctionalLink l,next;
+
+  PetscFunctionBeginUser;
+  if (!link) PetscFunctionReturn(0);
+  l     = *link;
+  *link = NULL;
+  for (; l; l=next) {
+    next = l->next;
+    ierr = PetscFree(l->name);CHKERRQ(ierr);
+    ierr = PetscFree(l);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SetInitialCondition"
+PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
+{
+  PetscScalar       *x;
+  PetscInt          cStart, cEnd, cEndInterior = user->cEndInterior, c;
+  PetscErrorCode    ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = VecGetArray(X, &x);CHKERRQ(ierr);
+  for (c = cStart; c < cEndInterior; ++c) {
+    PetscScalar    *xc;
+    ierr = DMPlexPointGlobalRef(dm,c,x,&xc);CHKERRQ(ierr);
+    if (xc) *xc = 0.; /* process ghost */
+  }
+  ierr = VecRestoreArray(X, &x);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "OutputVTK"
+static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);CHKERRQ(ierr);
+  ierr = PetscViewerSetType(*viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
+  ierr = PetscViewerFileSetName(*viewer, filename);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MonitorVTK"
+static PetscErrorCode MonitorVTK(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx)
+{
+  User           user = (User)ctx;
+  DM             dm;
+  PetscViewer    viewer;
+  char           filename[PETSC_MAX_PATH_LEN],*ftable = NULL;
+  PetscReal      xnorm;
+  Vec            E=0;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = PetscObjectSetName((PetscObject) X, "solution");CHKERRQ(ierr);
+  ierr = VecGetDM(X,&dm);CHKERRQ(ierr);
+  ierr = VecNorm(X,NORM_INFINITY,&xnorm);CHKERRQ(ierr);
+  ierr = VecDuplicate(X,&E);CHKERRQ(ierr);
+  if (stepnum >= 0) {           /* No summary for final time */
+    Model             mod = user->model;
+    PetscInt          c,cStart,cEnd,fcount,i;
+    size_t            ftableused,ftablealloc;
+    const PetscScalar *cellgeom_v,*x_v;
+    PetscScalar       *func_v;
+    DM                dmCell;
+    PetscReal         *fmin,*fmax,*fintegral,*ftmp;
+    fcount = mod->maxComputed+1;
+    ierr   = PetscMalloc4(fcount,PetscReal,&fmin,fcount,PetscReal,&fmax,fcount,PetscReal,&fintegral,fcount,PetscReal,&ftmp);CHKERRQ(ierr);
+    for (i=0; i<fcount; i++) {
+      fmin[i]      = PETSC_MAX_REAL;
+      fmax[i]      = PETSC_MIN_REAL;
+      fintegral[i] = 0;
+    }
+    ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
+    ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);
+    ierr = VecGetArrayRead(user->cellgeom,&cellgeom_v);CHKERRQ(ierr);
+    ierr = VecGetArrayRead(X,&x_v);CHKERRQ(ierr);
+    ierr = VecGetArray(E,&func_v);CHKERRQ(ierr);
+    for (c=cStart; c<user->cEndInterior; c++) {
+      const CellGeom    *cg;
+      const PetscScalar *cx;
+      PetscScalar *cfunc;
+      ierr = DMPlexPointLocalRead(dmCell,c,cellgeom_v,&cg);CHKERRQ(ierr);
+      ierr = DMPlexPointGlobalRead(dm,c,x_v,&cx);CHKERRQ(ierr);
+      ierr = DMPlexPointGlobalRef(dm,c,func_v,&cfunc);CHKERRQ(ierr);
+      if (!cx) continue;        /* not a global cell */
+      for (i=0; i<mod->numCall; i++) {
+        FunctionalLink flink = mod->functionalCall[i];
+        ierr = (*flink->func)(mod,time,cg->centroid,cx,&ftmp[i],flink->ctx);CHKERRQ(ierr);
+        *cfunc = ftmp[i]; /* save last functional */
+        if (i==mod->numCall-1) ftmp[i] = PetscAbs(ftmp[i]); /* allready a real, this assumes we want an abs - error */
+      }
+      for (i=0; i<fcount; i++) {
+        fmin[i]       = PetscMin(fmin[i],ftmp[i]);
+        fmax[i]       = PetscMax(fmax[i],ftmp[i]);
+        fintegral[i] += cg->volume * ftmp[i];
+      }
+    }
+    ierr = VecRestoreArrayRead(user->cellgeom,&cellgeom_v);CHKERRQ(ierr);
+    ierr = VecRestoreArrayRead(X,&x_v);CHKERRQ(ierr);
+    ierr = VecRestoreArray(E,&func_v);CHKERRQ(ierr);
+    ierr = MPI_Allreduce(MPI_IN_PLACE,fmin,fcount,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
+    ierr = MPI_Allreduce(MPI_IN_PLACE,fmax,fcount,MPIU_REAL,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
+    ierr = MPI_Allreduce(MPI_IN_PLACE,fintegral,fcount,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
+    if (fcount) { user->errorInf[user->ilev] = fmax[mod->numCall-1]; user->error2[user->ilev] = fintegral[mod->numCall-1]; } /* keep for convergence test */
+    ftablealloc = fcount * 100;
+    ftableused  = 0;
+    ierr        = PetscMalloc(ftablealloc,&ftable);CHKERRQ(ierr);
+    for (i=0; i<mod->numMonitored; i++) {
+      size_t         countused;
+      char           buffer[256],*p;
+      FunctionalLink flink = mod->functionalMonitored[i];
+      PetscInt       id    = flink->offset;
+      if (i % 3) {
+        ierr = PetscMemcpy(buffer,"  ",2);CHKERRQ(ierr);
+        p    = buffer + 2;
+      } else if (i) {
+        char newline[] = "\n";
+        ierr = PetscMemcpy(buffer,newline,sizeof newline-1);CHKERRQ(ierr);
+        p    = buffer + sizeof newline - 1;
+      } else {
+        p = buffer;
+      }
+      ierr = PetscSNPrintfCount(p,sizeof buffer-(p-buffer),"%12s functional: [%10.7G,%10.7G] integral=%10.7G",&countused,flink->name,fmin[id],fmax[id],fintegral[id]);CHKERRQ(ierr);
+      countused += p - buffer;
+      if (countused > ftablealloc-ftableused-1) { /* reallocate */
+        char *ftablenew;
+        ftablealloc = 2*ftablealloc + countused;
+        ierr = PetscMalloc(ftablealloc,&ftablenew);CHKERRQ(ierr);
+        ierr = PetscMemcpy(ftablenew,ftable,ftableused);CHKERRQ(ierr);
+        ierr = PetscFree(ftable);CHKERRQ(ierr);
+        ftable = ftablenew;
+      }
+      ierr = PetscMemcpy(ftable+ftableused,buffer,countused);CHKERRQ(ierr);
+      ftableused += countused;
+      ftable[ftableused] = 0;
+    }
+    ierr = PetscFree4(fmin,fmax,fintegral,ftmp);CHKERRQ(ierr);
+
+    ierr = PetscPrintf(PetscObjectComm((PetscObject)ts),"% 3D  time %8.4G  |x| %8.4G  %s\n",stepnum,time,xnorm,ftable ? ftable : "");CHKERRQ(ierr);
+    ierr = PetscFree(ftable);CHKERRQ(ierr);
+  }
+  if (user->vtkInterval < 1) {
+    ierr = VecDestroy(&E);CHKERRQ(ierr);
+    PetscFunctionReturn(0);
+  }
+  if (user->model->functionalCall && ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0))) {
+    Model          mod   = user->model;
+    FunctionalLink flink = mod->functionalCall[mod->numCall-1];
+    if (stepnum == -1) {        /* Final time is not multiple of normal time interval, write it anyway */
+      ierr = TSGetTimeStepNumber(ts,&stepnum);CHKERRQ(ierr);
+    }
+    ierr = PetscSNPrintf(filename,sizeof filename,"ex32-%03D-u.vtu",stepnum);CHKERRQ(ierr);
+    ierr = OutputVTK(dm,filename,&viewer);CHKERRQ(ierr);
+    ierr = VecView(X,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject) E, flink->name);CHKERRQ(ierr);
+    ierr = PetscSNPrintf(filename,sizeof filename,"ex32-%03D-%s.vtu",stepnum,flink->name);CHKERRQ(ierr);
+    ierr = OutputVTK(dm,filename,&viewer);CHKERRQ(ierr);
+    ierr = VecView(E,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+  }
+  ierr = VecDestroy(&E);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CreatePartitionVec"
+PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
+{
+  PetscSF        sfPoint;
+  PetscSection   coordSection;
+  Vec            coordinates;
+  PetscSection   sectionCell;
+  PetscScalar    *part;
+  PetscInt       cStart, cEnd, c;
+  PetscMPIInt    rank;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+  ierr = DMClone(dm, dmCell);CHKERRQ(ierr);
+  ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr);
+  ierr = DMSetPointSF(*dmCell, sfPoint);CHKERRQ(ierr);
+  ierr = DMPlexSetCoordinateSection(*dmCell, coordSection);CHKERRQ(ierr);
+  ierr = DMSetCoordinatesLocal(*dmCell, coordinates);CHKERRQ(ierr);
+  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    ierr = PetscSectionSetDof(sectionCell, c, 1);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(*dmCell, sectionCell);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&sectionCell);CHKERRQ(ierr);
+  ierr = DMCreateLocalVector(*dmCell, partition);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject)*partition, "partition");CHKERRQ(ierr);
+  ierr = VecGetArray(*partition, &part);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    PetscScalar *p;
+
+    ierr = DMPlexPointLocalRef(*dmCell, c, part, &p);CHKERRQ(ierr);
+    p[0] = rank;
+  }
+  ierr = VecRestoreArray(*partition, &part);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ApplyBC"
+static PetscErrorCode ApplyBC(DM dm, PetscReal time, Vec locX, User user)
+{
+  Model             mod = user->model;
+  const char        *name = "Face Sets";
+  DM                dmFace;
+  IS                idIS;
+  const PetscInt    *ids;
+  PetscScalar       *x;
+  const PetscScalar *facegeom;
+  PetscInt          numFS, fs;
+  PetscErrorCode    ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);
+  ierr = DMPlexGetLabelIdIS(dm, name, &idIS);CHKERRQ(ierr);
+  if (!idIS) PetscFunctionReturn(0);
+  ierr = ISGetLocalSize(idIS, &numFS);CHKERRQ(ierr);
+  ierr = ISGetIndices(idIS, &ids);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(user->facegeom, &facegeom);CHKERRQ(ierr);
+  ierr = VecGetArray(locX, &x);CHKERRQ(ierr);
+  for (fs = 0; fs < numFS; ++fs) {
+    BoundaryFunction bcFunc;
+    void             *bcCtx;
+    IS               faceIS;
+    const PetscInt   *faces;
+    PetscInt         numFaces, f;
+    ierr = ModelBoundaryFind(mod,ids[fs],&bcFunc,&bcCtx);CHKERRQ(ierr);
+    ierr = DMPlexGetStratumIS(dm, name, ids[fs], &faceIS);CHKERRQ(ierr);
+    ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
+    ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
+    for (f = 0; f < numFaces; ++f) {
+      const PetscInt    face = faces[f], *cells;
+      const PetscScalar *xI;
+      PetscScalar       *xG;
+      const FaceGeom    *fg;
+      ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr);
+      ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
+      ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr);
+      ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr);
+      ierr = (*bcFunc)(mod, time, fg->centroid, fg->normal, xI, xG, bcCtx);CHKERRQ(ierr);
+    }
+    ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr);
+    ierr = ISDestroy(&faceIS);CHKERRQ(ierr);
+  }
+  ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr);
+  ierr = VecRestoreArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(idIS, &ids);CHKERRQ(ierr);
+  ierr = ISDestroy(&idIS);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ApplyLaplacianLocal"
+static PetscErrorCode ApplyLaplacianLocal(DM dm,DM dmFace,DM dmCell,PetscReal time,Vec locX,Vec F,User user)
+{
+  /* Physics           phys = user->model->physics; */
+  DMLabel           ghostLabel;
+  PetscErrorCode    ierr;
+  const PetscScalar *facegeom, *cellgeom, *x;
+  PetscScalar       *f;
+  PetscInt          fStart, fEnd, face;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(locX,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
+  ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
+  for (face = fStart; face < fEnd; ++face) {
+    const PetscInt    *cells;
+    PetscInt          ghost;
+    PetscScalar       *fL,*fR,v[DIM],dx,area,flux,alpha;
+    const FaceGeom    *fg;
+    const CellGeom    *cgL,*cgR;
+    const PetscScalar *xL,*xR;
+
+    ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr);
+    if (ghost >= 0) continue;
+    ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
+    ierr = DMPlexPointLocalRead(dmFace,face,facegeom,&fg);CHKERRQ(ierr);
+    ierr = DMPlexPointLocalRead(dmCell,cells[0],cellgeom,&cgL);CHKERRQ(ierr);
+    ierr = DMPlexPointLocalRead(dmCell,cells[1],cellgeom,&cgR);CHKERRQ(ierr);
+    ierr = DMPlexPointLocalRead(dm,cells[0],x,&xL);CHKERRQ(ierr);
+    ierr = DMPlexPointLocalRead(dm,cells[1],x,&xR);CHKERRQ(ierr);
+    ierr = DMPlexPointGlobalRef(dm,cells[0],f,&fL);CHKERRQ(ierr);
+    ierr = DMPlexPointGlobalRef(dm,cells[1],f,&fR);CHKERRQ(ierr);
+    ierr = (*user->model->alpha)(user->model,time,fg->centroid,&alpha,user->model->physics);CHKERRQ(ierr); /* store this? */
+    WaxpyD(-1,cgL->centroid,cgR->centroid,v); /* dx vector */
+    dx = Norm2D(v);
+    area = Norm2D(fg->normal);
+    flux = alpha*(xR[0]-xL[0])/dx;
+    if (fL) fL[0] -= area*flux / cgL->volume;
+    if (fR) fR[0] += area*flux / cgR->volume;
+  }
+  ierr = VecRestoreArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
+  ierr = VecRestoreArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);
+  ierr = VecRestoreArrayRead(locX,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ComputeFunctionLap"
+/*
+ Finite difference operator for Jacobian.  The operator in SNES form.
+ */
+PetscErrorCode ComputeFunctionLap(SNES snes,Vec U,Vec FU,void *ctx)
+{
+  User           user = (User)ctx;
+  DM             dm, dmFace, dmCell;
+  Vec            locU;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
+  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);
+  ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);
+  ierr = DMGetLocalVector(dm,&locU);CHKERRQ(ierr);
+  ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, locU);CHKERRQ(ierr); /* get ghost values - echange */
+  ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, locU);CHKERRQ(ierr);
+
+  ierr = ApplyBC(dm, 0., locU, user);CHKERRQ(ierr); /* set ghost values - boundary */
+
+  ierr = VecZeroEntries(FU);CHKERRQ(ierr);
+  ierr = ApplyLaplacianLocal(dm,dmFace,dmCell,0.,locU,FU,user);CHKERRQ(ierr);
+
+  ierr = DMRestoreLocalVector(dm,&locU);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "RHSFunctionLap"
+/*
+ Laplacian + Source term for TS RHS: Au - b
+ */
+PetscErrorCode RHSFunctionLap(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
+{
+  PetscErrorCode ierr;
+  Mat            A;
+  User           user = (User)ctx;
+  Model          mod = user->model;
+  PetscInt       cStart, cEnd, c;
+  DM             dm,dmCell;
+  const PetscScalar *cgeom_v;
+  PetscScalar *x_v,source[1];
+
+  PetscFunctionBeginUser;
+  ierr = TSGetRHSJacobian(ts,&A,NULL,NULL,&ctx);CHKERRQ(ierr);
+  ierr = MatMult(A,globalin,globalout);CHKERRQ(ierr);
+  /* subtract off RHS to solve Ax=b */
+  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
+  ierr = VecGetDM(user->cellgeom, &dmCell);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(user->cellgeom, &cgeom_v);CHKERRQ(ierr);
+  ierr = VecGetArray(globalout, &x_v);CHKERRQ(ierr);
+  for (c = cStart; c < user->cEndInterior; ++c) {
+    CellGeom *cg;
+    PetscScalar *xc;
+    ierr = DMPlexPointLocalRead(dmCell,c,cgeom_v,&cg);CHKERRQ(ierr);
+    ierr = (*mod->source)(mod,t,cg->centroid,source,mod->physics);CHKERRQ(ierr);
+    ierr = DMPlexPointGlobalRef(dm,c,x_v,&xc);CHKERRQ(ierr);
+    if (xc) *xc -= source[0];
+  }
+  ierr = VecRestoreArrayRead(user->cellgeom, &cgeom_v);CHKERRQ(ierr);
+  ierr = VecRestoreArray(globalout, &x_v);CHKERRQ(ierr);
+  VecScale(globalout, -1.); /* why is this needed to get SNES residual correct? */
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char **argv)
+{
+  MPI_Comm          comm;
+  User              user;
+  Model             mod;
+  Physics           phys;
+  DM                dm, dmDist;
+  PetscReal         ftime,dt;
+  PetscInt          dim, overlap, nsteps, nrefine,i,ilev;
+  int               CPU_word_size = 0, IO_word_size = 0, exoid;
+  float             version;
+  TS                ts;
+  TSConvergedReason reason;
+  Vec               X;
+  Mat               A;
+  PetscViewer       viewer;
+  PetscMPIInt       rank;
+  char              filename[PETSC_MAX_PATH_LEN] = "./square-4x4.exo";
+  PetscBool         vtkCellGeom;
+  PetscErrorCode    ierr;
+
+  ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);
+  comm = PETSC_COMM_WORLD;
+  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
+
+  ierr = PetscNew(struct _n_User,&user);CHKERRQ(ierr);
+  ierr = PetscNew(struct _n_Model,&user->model);CHKERRQ(ierr);
+  ierr = PetscNew(struct _n_Physics,&user->model->physics);CHKERRQ(ierr);
+  mod  = user->model;
+  phys = mod->physics;
+  mod->comm = comm;
+
+  ierr = PetscOptionsBegin(comm,NULL,"Unstructured Finite Volume Options","");CHKERRQ(ierr);
+  {
+    ierr              = PetscOptionsString("-f","Exodus.II filename to read","",filename,filename,sizeof(filename),NULL);CHKERRQ(ierr);
+    user->vtkInterval = 1;
+    ierr = PetscOptionsInt("-ufv_vtk_interval","VTK output interval (0 to disable)","",user->vtkInterval,&user->vtkInterval,NULL);CHKERRQ(ierr);
+    overlap = 1;
+    ierr = PetscOptionsInt("-ufv_mesh_overlap","Number of cells to overlap partitions","",overlap,&overlap,NULL);CHKERRQ(ierr);
+    vtkCellGeom = PETSC_FALSE;
+    nrefine=0;
+    ierr = PetscOptionsInt("-nrefine","Number of refinment steps","",nrefine,&nrefine,NULL);CHKERRQ(ierr);
+    if (nrefine>=32) nrefine = 31;
+    ierr = PetscOptionsBool("-ufv_vtk_cellgeom","Write cell geometry (for debugging)","",vtkCellGeom,&vtkCellGeom,NULL);CHKERRQ(ierr);
+    ierr = PetscMemzero(phys,sizeof(struct _n_Physics));CHKERRQ(ierr);
+    ierr = PhysicsCreate_Lap(mod,phys);CHKERRQ(ierr);
+    /* Count number of fields and dofs */
+    for (phys->nfields=0,phys->dof=0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
+
+    if (phys->dof != 1) SETERRQ1(comm,PETSC_ERR_ARG_WRONGSTATE,"Physics dof=%d",phys->dof);
+    ierr = ModelFunctionalSetFromOptions(mod);CHKERRQ(ierr);
+  }
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  ierr = SetUpBoundaries(comm, user);CHKERRQ(ierr);
+
+  /* outer congerance loop */
+  for (ilev=0;ilev<=nrefine;ilev++) {
+    user->ilev = ilev;
+    /* create base mesh */
+    if (!rank) {
+      exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
+      if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ex_open(\"%s\",...) did not return a valid file ID",filename);
+    } else exoid = -1;                 /* Not used */
+    ierr = DMPlexCreateExodus(comm, exoid, PETSC_TRUE, &dm);CHKERRQ(ierr);
+    if (!rank) {ierr = ex_close(exoid);CHKERRQ(ierr);}
+
+    /* refine mesh */
+    /* ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); */
+    for (i=0;i<ilev;i++) {
+      DM refinedMesh     = NULL;
+      /* Refine mesh using ... */
+      ierr = DMRefine(dm, comm, &refinedMesh);CHKERRQ(ierr);
+      if (refinedMesh) {
+        DMLabel  faceSets;
+        PetscInt fStart, fEnd;
+        ierr = DMDestroy(&dm);CHKERRQ(ierr);
+        ierr = DMPlexGetHeightStratum(refinedMesh, 1, &fStart, &fEnd);CHKERRQ(ierr);
+        ierr = DMPlexGetLabel(refinedMesh, "Face Sets", &faceSets);CHKERRQ(ierr);
+        if (faceSets) {ierr = DMLabelFilter(faceSets, fStart, fEnd);CHKERRQ(ierr);}
+        dm  = refinedMesh;
+      }
+    }
+
+    /* Distribute mesh */
+    ierr = DMPlexDistribute(dm, "metis", overlap, NULL, &dmDist);CHKERRQ(ierr);
+    if (dmDist) {
+      ierr = DMDestroy(&dm);CHKERRQ(ierr);
+      dm   = dmDist;
+    }
+    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
+
+    {
+      DM gdm;
+      ierr = DMPlexGetHeightStratum(dm, 0, NULL, &user->cEndInterior);CHKERRQ(ierr);
+      ierr = DMPlexConstructGhostCells(dm, NULL, &user->numGhostCells, &gdm);CHKERRQ(ierr);
+      ierr = DMDestroy(&dm);CHKERRQ(ierr);
+      dm   = gdm;
+    }
+    /* DMView(dm,PETSC_VIEWER_STDOUT_WORLD); */
+    ierr = ConstructGeometry(dm, &user->facegeom, &user->cellgeom, user);CHKERRQ(ierr);
+    if (0) {ierr = VecView(user->cellgeom, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
+    ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+    ierr = DMPlexSetPreallocationCenterDimension(dm, 0);CHKERRQ(ierr);
+
+    /* Set up DM with section describing local vector and configure local vector. */
+    ierr = SetUpLocalSpace(dm, user);CHKERRQ(ierr);
+
+    ierr = DMCreateGlobalVector(dm, &X);CHKERRQ(ierr);
+    ierr = SetInitialCondition(dm, X, user);CHKERRQ(ierr);
+
+    if (vtkCellGeom) {
+      DM  dmCell;
+      Vec partition;
+
+      ierr = OutputVTK(dm, "ex32-cellgeom.vtk", &viewer);CHKERRQ(ierr);
+      ierr = VecView(user->cellgeom, viewer);CHKERRQ(ierr);
+      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+      ierr = CreatePartitionVec(dm, &dmCell, &partition);CHKERRQ(ierr);
+      ierr = OutputVTK(dmCell, "ex32-partition.vtk", &viewer);CHKERRQ(ierr);
+      ierr = VecView(partition, viewer);CHKERRQ(ierr);
+      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+      ierr = VecDestroy(&partition);CHKERRQ(ierr);
+      ierr = DMDestroy(&dmCell);CHKERRQ(ierr);
+    }
+
+    ierr = DMCreateMatrix(dm, &A);CHKERRQ(ierr);
+    { /* finite difference Jacobian A */
+      SNES snes;
+      MatStructure flag;
+
+      ierr = SNESCreate(comm, &snes);CHKERRQ(ierr);
+      ierr = SNESSetDM(snes,dm);CHKERRQ(ierr);
+      ierr = SNESSetFunction(snes,NULL,ComputeFunctionLap,user);CHKERRQ(ierr);
+      ierr = SNESSetJacobian(snes,A,A,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);
+
+      ierr = SNESComputeJacobian(snes,X,&A,&A,&flag);CHKERRQ(ierr);
+      /* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
+      ierr = SNESDestroy(&snes);CHKERRQ(ierr);
+    }
+    ierr = MatGetSize(A, &user->N[ilev], NULL);CHKERRQ(ierr);
+
+    /* create TS */
+    ierr = TSCreate(comm, &ts);CHKERRQ(ierr);
+    ierr = TSSetProblemType(ts,TS_LINEAR);CHKERRQ(ierr);
+    ierr = TSSetType(ts, TSBEULER);CHKERRQ(ierr);
+    ierr = TSSetDM(ts, dm);CHKERRQ(ierr);
+    ierr = TSMonitorSet(ts,MonitorVTK,user,NULL);CHKERRQ(ierr);
+    ierr = TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&user);CHKERRQ(ierr);
+    ierr = TSSetRHSFunction(ts,NULL,RHSFunctionLap,user);CHKERRQ(ierr);
+    ierr = TSSetDuration(ts,1,2.0);CHKERRQ(ierr);
+    dt   = 1.e20;
+    ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
+    ierr = TSSetSolution(ts,X);CHKERRQ(ierr);
+    ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+    /* solve */
+    ierr = TSSolve(ts,X);CHKERRQ(ierr);
+    ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
+    ierr = TSGetTimeStepNumber(ts,&nsteps);CHKERRQ(ierr);
+    ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);
+    PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps. Errors: |e|_inf=%e |e|_2=%e, N=%d\n",TSConvergedReasons[reason],ftime,nsteps,user->errorInf[ilev],user->error2[ilev],user->N[ilev]);
+
+    /* clean up */
+    ierr = TSDestroy(&ts);CHKERRQ(ierr);
+    ierr = DMDestroy(&dm);CHKERRQ(ierr);
+    ierr = VecDestroy(&user->cellgeom);CHKERRQ(ierr);
+    ierr = VecDestroy(&user->facegeom);CHKERRQ(ierr);
+    ierr = VecDestroy(&X);CHKERRQ(ierr);
+    ierr = MatDestroy(&A);CHKERRQ(ierr);
+  }
+
+  /* print out convergence rates */
+  if (nrefine>0 && mod->numCall>0) {
+    PetscPrintf(PETSC_COMM_WORLD,"\tN | Convergance rates |e|_inf\t\t|e|_2  (%d)\n",mod->numCall);
+    for (ilev=1;ilev<=nrefine;ilev++) {
+      const PetscReal log2r = 1.0/log(2.0);
+      PetscReal ratio = user->errorInf[ilev-1]/user->errorInf[ilev];
+      const PetscReal convergeRate_Inf = log(ratio)*log2r;
+      ratio = user->error2[ilev-1]/user->error2[ilev];
+      const PetscReal convergeRate_2 = log(ratio)*log2r;
+      PetscPrintf(PETSC_COMM_WORLD,"\t%d\t\t\t%e\t%e\n",user->N[ilev],convergeRate_Inf,convergeRate_2);
+    }
+  }
+
+  ierr = FunctionalLinkDestroy(&user->model->functionalRegistry);CHKERRQ(ierr);
+  ierr = BoundaryLinkDestroy(&user->model->boundary);CHKERRQ(ierr);
+  ierr = PetscFree(user->model->functionalMonitored);CHKERRQ(ierr);
+  ierr = PetscFree(user->model->functionalCall);CHKERRQ(ierr);
+  ierr = PetscFree(user->model->physics->data);CHKERRQ(ierr);
+  ierr = PetscFree(user->model->physics);CHKERRQ(ierr);
+  ierr = PetscFree(user->model);CHKERRQ(ierr);
+  ierr = PetscFree(user);CHKERRQ(ierr);
+
+  /* { */
+  /*   // setup interpolation */
+  /*   const PetscInt nref=3; */
+  /*   PetscReal targetx[nref+1],sourcex[nref+1],RCoef[nref*nref]; */
+  /*   int i,j,k; */
+  /*   for (i=0,j=-(nref-1)/2;i<=nref;i++,j++) targetx[i] = j; */
+  /*   //  normal interpolant */
+  /*   //  sources:     0             1            2           3 */
+  /*   //  targets (*)  |      o      |      *     |     o     | */
+  /*   //  coord       -1             0.           1           2 */
+  /*   for (i=0,j=-(nref-1)/2;i<=nref;i++,j++) sourcex[i] = j; */
+  /*   printf("source: "); */
+  /*   for (i=0;i<=nref;i++) printf(" %e",sourcex[i]); printf("\ntargetx:"); */
+  /*   for (i=0;i<=nref;i++) printf(" %e",targetx[i]); printf("\n"); */
+  /*   // PetscDTReconstructPoly(degree,nsource,*sourcex,ntarget,*targetx,*R) */
+  /*   ierr = PetscDTReconstructPoly(nref-1,nref,sourcex,nref,targetx,RCoef);CHKERRQ(ierr); */
+  /*   for (i=0,k=0;i<nref;i++) { */
+  /*     printf("RCoef:"); */
+  /*     for (j=0;j<nref;j++,k++) printf(" %e",RCoef[k]); */
+  /*     printf("\n"); */
+  /*   } */
+  /* } */
+
+  ierr = PetscFinalize();
+  return(0);
+  }

# File src/ts/examples/tutorials/makefile

 	-${CLINKER} -o ex31 ex31.o${PETSC_TS_LIB}
 	${RM} ex31.o   +ex33: ex33.o chkopts + -${CLINKER} -o ex33 ex33.o ${PETSC_TS_LIB} +${RM} ex33.o
+
 #---------------------------------------------------------------------------------
 runex1:
 	-@${MPIEXEC} -n 1 ./ex1 -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_monitor_pseudo > ex1_1.tmp 2>&1; \ ${DIFF} output/ex31_1.out ex31_1.tmp || echo  ${PWD} "\nPossible problem with ex20_1, diffs above \n========================================="; \ ${RM} -f ex31_1.tmp
 
+runex33:
+	-@${MPIEXEC} -n 1 ./ex33 -nrefine 3 -f${PETSC_DIR}/share/petsc/datafiles/meshes/square-4x4.exo -monitor error > ex33_1.tmp 2>&1; \
+	   ${DIFF} output/ex33_1.out ex33_1.tmp || echo${PWD} "\nPossible problem with ex20_1, diffs above \n=========================================";  \
+	   \${RM} -f ex33_1.tmp
+
 TESTEXAMPLES_C		  = ex1.PETSc runex1 runex1_2 ex1.rm ex2.PETSc runex2 ex2.rm ex3.PETSc runex3 runex3_2 ex3.rm \
                             ex4.PETSc runex4 runex4_2 runex4_3 runex4_4 ex4.rm ex5.PETSc ex5.rm \
                             ex6.PETSc runex6 ex6.rm ex7.PETSc runex7 runex7_2 runex7_3 ex7.rm \

# File src/ts/examples/tutorials/output/ex33_1.out

+  0  time        0  |x|        0         error functional: [0.0002365708,0.05666167] integral=0.01916224
+  1  time    1e+20  |x|  0.06243         error functional: [8.420959e-05,0.01281517] integral=0.004089645
+CONVERGED_ITS at time 1e+20 after 1 steps. Errors: |e|_inf=1.281517e-02 |e|_2=4.089645e-03, N=16
+  0  time        0  |x|        0         error functional: [1.513981e-05,0.06212672] integral=0.01812478
+  1  time    1e+20  |x|  0.06421         error functional: [7.53915e-05,0.003876525] integral=0.001104024
+CONVERGED_ITS at time 1e+20 after 1 steps. Errors: |e|_inf=3.876525e-03 |e|_2=1.104024e-03, N=64
+  0  time        0  |x|        0         error functional: [9.518126e-07,0.06236227] integral=0.01786457
+  1  time    1e+20  |x|  0.06296         error functional: [8.375495e-07,0.001091853] integral=0.0002809029
+CONVERGED_ITS at time 1e+20 after 1 steps. Errors: |e|_inf=1.091853e-03 |e|_2=2.809029e-04, N=256
+	N | Convergance rates |e|_inf		|e|_2  (1)
+	64			1.725016e+00	1.889204e+00
+	256			1.827986e+00	1.974628e+00`