# Commits

committed 73c2058

TSSYMPEULER and examples

Adds a new TS implementation, TSSYMPEULER, which implements
a semi-implicit symplectic integration scheme. In its current
form it is mostly useful as a demonstration of using a partitioned
righthand side with a structured integrator. Two new TS tutorials
are included.

• Participants
• Parent commits 8294b37
• Branches psanan/ts-rhs-partitions

# File include/petscts.h

• Ignore whitespace
 #define TSARKIMEX         "arkimex"
 #define TSROSW            "rosw"
 #define TSEIMEX           "eimex"
+#define TSSYMPEULER       "sympeuler"
 /*E
     TSProblemType - Determines the type of problem this TS object is to be used to solve


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

• Ignore whitespace
+static const char help[] = "Demonstrate the interface to solve DAE systems with partitioned RHS functions\n";
+/*
+   Consider the following second order ODE:
+
+   q''(t) + q(t) = \epsilon f(q)
+
+   where \epsilon may or may not be small
+
+   As a first order system, we  have
+
+   q'(t) = p(t)
+   p'(t) = -q(t) + \epsilon f(q)
+
+   (We choose f(q) = -q^3)
+
+   Write X = [p; q] and partition the righthand side
+
+   dX/dt =(d/dt) [q; p] =  [p(t); 0] + [0; -q(t)] + [0 ; \epsilon f(q)]
+
+                           -- P ---   -------------- Q ----------------
+                           ------- Fast ---------   ------ Slow -------
+
+  Different integration schemes would split the same RHS differently,
+  perhaps based on the value of \epsilon. Here, we use the first splitting
+  with a simple symplectic integrator. A future exponential integrator might
+  use the second splitting.
+
+  Accepts an extra option -epsilon
+*/
+
+#include <petscts.h>
+
+/* User-defined Functions */
+static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
+static PetscErrorCode FormRHSFunctionP(TS,PetscReal,Vec,Vec,void*);
+static PetscErrorCode FormRHSFunctionQ(TS,PetscReal,Vec,Vec,void*);
+static PetscErrorCode FormJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
+static PetscErrorCode FormJacobianP(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
+static PetscErrorCode FormJacobianQ(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
+static PetscErrorCode FormInitialSolution(TS,Vec,void*);
+
+/* User context */
+typedef struct
+{
+  PetscReal epsilon;
+} User;
+
+/*
+Main function
+*/
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char* argv[])
+{
+  PetscErrorCode    ierr;
+  PetscMPIInt       size;
+  Vec               X;
+  Mat               J;
+  TS                ts;
+  DM                dm;
+  User              user;
+  PetscScalar       epsilon = 1;
+  PetscReal         dt = 0.1,maxtime = 50,ftime;
+  PetscInt          steps;
+  TSConvergedReason reason;
+
+  ierr = PetscInitialize(&argc, &argv, (char*) 0,help);CHKERRQ(ierr);
+  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);
+  if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only");
+
+  ierr = PetscOptionsGetScalar(NULL,"-epsilon",&epsilon,NULL);CHKERRQ(ierr);
+
+  /* User Parameters */
+  user.epsilon = epsilon;
+
+  /* Create TS */
+  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
+  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
+
+  /* Create Solution Vector and Jacobian Matrix*/
+  ierr = VecCreateSeq(PETSC_COMM_WORLD,2,&X);CHKERRQ(ierr);
+  ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&J);CHKERRQ(ierr);
+
+  /* Set Initial Conditions and Parameters*/
+  ierr = FormInitialSolution(ts,X,&user);
+  ierr = TSSetSolution(ts,X);CHKERRQ(ierr);
+  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,PETSC_MAX_INT,maxtime);
+  ierr = TSSetType(ts,TSSYMPEULER);CHKERRQ(ierr);
+
+  /* Register Partitioned RHS Functions and Jacobians */
+  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr);
+  ierr = DMTSSetRHSPartitionFunction(dm,TS_SYMP_PARTITION,TS_SYMP_P_SLOT,FormRHSFunctionP,&user);CHKERRQ(ierr);
+  ierr = DMTSSetRHSPartitionFunction(dm,TS_SYMP_PARTITION,TS_SYMP_Q_SLOT,FormRHSFunctionQ,&user);CHKERRQ(ierr);
+  ierr = TSSetRHSJacobian(ts,J,J,FormJacobian,&user);CHKERRQ(ierr);
+  ierr = DMTSSetRHSPartitionJacobian(dm,TS_SYMP_PARTITION,TS_SYMP_Q_SLOT,FormJacobianQ,&user);CHKERRQ(ierr);
+  ierr = DMTSSetRHSPartitionJacobian(dm,TS_SYMP_PARTITION,TS_SYMP_P_SLOT,FormJacobianP,&user);CHKERRQ(ierr);
+
+  /* Set from options */
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+  /* Solve */
+  ierr = TSSolve(ts,X); CHKERRQ(ierr);
+  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
+  ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);
+  ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,steps);CHKERRQ(ierr);
+
+  /* Clean Up */
+  ierr = MatDestroy(&J);CHKERRQ(ierr);
+  ierr = VecDestroy(&X);CHKERRQ(ierr);
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+
+  return EXIT_SUCCESS;
+}
+
+/*
+Initial Conditions
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormInitialSolution"
+static PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
+{
+  PetscScalar    *x;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  x[0] = 1;
+  x[1] = 1;
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS function (complete)
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormRHSFunction"
+static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void* ctx)
+{
+  User             *user = (User*) ctx;
+  PetscScalar      *x,*f;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
+  f[0] = x[1];
+  f[1] = -x[0] -(user->epsilon*x[0]*x[0]*x[0]);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS function for 'P' terms
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormRHSFunctionP"
+static PetscErrorCode FormRHSFunctionP(TS ts, PetscReal t, Vec X, Vec F, void* ctx)
+{
+  User            *user = (User*) ctx;
+  PetscScalar     *x,*f;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
+  f[0] = 0;
+  f[1] = -x[0] -(user->epsilon*x[0]*x[0]*x[0]);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS function for 'Q' terms
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormRHSFunctionQ"
+static PetscErrorCode FormRHSFunctionQ(TS ts, PetscReal t, Vec X, Vec F, void* ctx)
+{
+  PetscScalar     *x,*f;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
+  f[0] = x[1];
+  f[1] = 0;
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS Jacobian (full)
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormJacobian"
+PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ctx)
+{
+  User              *user = (User*)ctx;
+  PetscErrorCode    ierr;
+  PetscScalar       v[4],*x;
+  PetscInt          idxm[2] = {0,1};
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  v[0] = 0;                                   v[1] = 1;
+  v[2] = -1 -(user->epsilon * x[0] * x[0]);   v[3] = 0;
+  ierr = MatSetValues(*B,2,idxm,2,idxm,v,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  if (*J != *B) {
+    ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+    ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  }
+  *flag = SAME_NONZERO_PATTERN;
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS Jacobian (P)
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormJacobianP"
+PetscErrorCode FormJacobianP(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ctx)
+{
+  User              *user = (User*)ctx;
+  PetscErrorCode    ierr;
+  PetscScalar       v[4], *x;
+  PetscInt          idxm[2] = {0,1};
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  v[0] = 0;                                   v[1] = 0;
+  v[2] = -1 -(user->epsilon * x[0] * x[0]);   v[3] = 0;
+  ierr = MatSetValues(*B,2,idxm,2,idxm,v,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  if (*J != *B) {
+    ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+    ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  }
+  *flag = SAME_NONZERO_PATTERN;
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS Jacobian (Q)
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormJacobianQ"
+PetscErrorCode FormJacobianQ(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ctx)
+{
+  PetscErrorCode    ierr;
+  PetscScalar       v[4], *x;
+  PetscInt          idxm[2] = {0,1};
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  v[0] = 0;   v[1] = 1;
+  v[2] = 0;   v[3] = 0;
+  ierr = MatSetValues(*B,2,idxm,2,idxm,v,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  if (*J != *B) {
+    ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+    ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  }
+  *flag = SAME_NONZERO_PATTERN;
+  PetscFunctionReturn(0);
+}

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

• Ignore whitespace
+static const char help[] = "Demonstrate the use of an integrator requiring a partitioned RHS\n";
+/*
+
+  Test for a Hamiltonian system with a variable mass matrix
+
+  This follows Hairer/Lubich/Wanner, "Geometric Numerical Integration", 2nd ed.,  ch XIII example 10.1
+
+  Coordinates are
+  q0 [angle]
+  q1 [r-1]
+  p0 = (1+q_1)^2 \dot q_0
+  p1 = \dot q-1
+
+  Note that this problem actually admits an explicit algorithm with an even finer splitting of the RHS (detailed in the book).
+  Here we use it as a test for the symplectic Euler method with a nontrivial explicit step.
+
+  Note that in the current implementation, the 'P' terms are treated implicitly and the 'Q' terms explicitly. Thus,
+  a Jacobian is only provided for one of the partitioned RHS functions.
+
+  Accepts an extra option -epsilon
+*/
+
+#include <petscts.h>
+
+/* User-defined Functions */
+static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
+static PetscErrorCode FormRHSFunctionP(TS,PetscReal,Vec,Vec,void*);
+static PetscErrorCode FormRHSFunctionQ(TS,PetscReal,Vec,Vec,void*);
+static PetscErrorCode  FormJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
+static PetscErrorCode  FormJacobianP(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
+static PetscErrorCode FormInitialSolution(TS,Vec,void*);
+
+/* User context */
+typedef struct
+{
+  PetscReal epsilon;
+} User;
+
+/*
+Main function
+*/
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char* argv[])
+{
+  PetscErrorCode      ierr;
+  PetscMPIInt         size;
+  Vec                 X;
+  Mat                 J;
+  TS                  ts;
+  DM                  dm;
+  User                user;
+  PetscScalar         epsilon = 0.1;
+  PetscReal           dt = 0.15,maxtime = 10,ftime;
+  PetscInt            steps;
+  TSConvergedReason   reason;
+
+  ierr = PetscInitialize(&argc, &argv, (char*) 0,help);CHKERRQ(ierr);
+  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);
+  if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only");
+
+  ierr = PetscOptionsGetScalar(NULL,"-epsilon",&epsilon,NULL);CHKERRQ(ierr);
+
+  /* User Parameters */
+  user.epsilon = epsilon;
+
+  /* Create TS */
+  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
+  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
+
+  /* Create Vector and Matrix for Solution and Jacobian */
+  ierr = VecCreateSeq(PETSC_COMM_WORLD,4,&X);CHKERRQ(ierr);
+  ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD,4,4,4,NULL,&J);CHKERRQ(ierr);
+
+  /* Set Initial Conditions and time Parameters*/
+  ierr = FormInitialSolution(ts,X,&user);
+  ierr = TSSetSolution(ts,X);CHKERRQ(ierr);
+  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,PETSC_MAX_INT,maxtime);
+
+  /* Set RHS Functions amd Jacobians*/
+  ierr = DMTSSetRHSFunction(dm,FormRHSFunction,&user);CHKERRQ(ierr);
+
+  /* Register Partitioned RHS Functions and Jacobians */
+  ierr = DMTSSetRHSPartitionFunction(dm,TS_SYMP_PARTITION,TS_SYMP_P_SLOT,FormRHSFunctionP,&user);CHKERRQ(ierr);
+  ierr = DMTSSetRHSPartitionFunction(dm,TS_SYMP_PARTITION,TS_SYMP_Q_SLOT,FormRHSFunctionQ,&user);CHKERRQ(ierr);
+  ierr = TSSetRHSJacobian(ts, J, J, FormJacobian, &user);CHKERRQ(ierr);
+  ierr = DMTSSetRHSPartitionJacobian(dm,TS_SYMP_PARTITION,TS_SYMP_P_SLOT, FormJacobianP, &user);CHKERRQ(ierr);
+
+  /* Set TS Type  */
+  ierr = TSSetType(ts,TSSYMPEULER);CHKERRQ(ierr);
+
+  /* Set from Options */
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+  /* Solve */
+  ierr = TSSolve(ts,X); CHKERRQ(ierr);
+  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
+  ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);
+  ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,steps);CHKERRQ(ierr);
+
+  /* Clean Up */
+  ierr = VecDestroy(&X);CHKERRQ(ierr);
+  ierr = MatDestroy(&J);CHKERRQ(ierr);
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+
+  return EXIT_SUCCESS;
+}
+
+/*
+Initial Conditions
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormInitialSolution"
+static PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
+{
+  PetscScalar    *x;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  x[0] = 1.57079632679;
+  x[1] = 0.2;
+  x[2] = 0.3;
+  x[3] = 0.1;
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS function (complete)
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormRHSFunction"
+static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void* ctx)
+{
+  User             *user = (User*) ctx;
+  PetscScalar      *x,*f,opq1;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
+  opq1 = 1 + x[1];
+  f[0] = x[2]/(opq1*opq1);
+  f[1] = x[3];
+  f[2] = -(1 + x[1]) * PetscSinScalar(x[0]);
+  f[3] = (x[2] * x[2])/(opq1 * opq1 * opq1) - x[1]/(user->epsilon * user->epsilon) + PetscCosScalar(x[0]);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS function for 'P' terms
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormRHSFunctionP"
+static PetscErrorCode FormRHSFunctionP(TS ts, PetscReal t, Vec X, Vec F, void* ctx)
+{
+  User             *user = (User*) ctx;
+  PetscScalar      *x,*f,opq1;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
+  opq1 = 1 + x[1];
+  f[0] = 0;
+  f[1] = 0;
+  f[2] = -(1 + x[1]) * PetscSinScalar(x[0]);
+  f[3] = (x[2] * x[2])/(opq1 * opq1 * opq1) - x[1]/(user->epsilon * user->epsilon) + PetscCosScalar(x[0]);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS function for 'Q' terms
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormRHSFunctionQ"
+static PetscErrorCode FormRHSFunctionQ(TS ts, PetscReal t, Vec X, Vec F, void* ctx)
+{
+  PetscScalar      *x,*f,opq1;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
+  opq1 = 1 + x[1];
+  f[0] = x[2]/(opq1*opq1);
+  f[1] = x[3];
+  f[2] = 0;
+  f[3] = 0;
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS Jacobian (full)
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormJacobian"
+PetscErrorCode FormJacobian(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ctx)
+{
+  User              *user = (User*)ctx;
+  PetscErrorCode    ierr;
+  PetscScalar       v[16],*x,opq1;
+  PetscInt          idxm[4] = {0,1,2,3};
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  opq1 = 1 + x[1];
+ v[0]  = 1/(opq1*opq1);           v[1]  = 0; v[2] = 0;                           v[3]  = -2/(opq1*opq1*opq1)*x[2];
+ v[4]  = 0;                       v[5]  = 1; v[6] = 0;                           v[5]  = 0;
+ v[8]  = 0;                       v[9]  = 0; v[10] = -opq1*PetscCosScalar(x[0]); v[9]  = -PetscSinScalar(x[0]);
+ v[12] = 2/(opq1*opq1*opq1)*x[2]; v[13] = 0; v[14] = -PetscSinScalar(x[0]);      v[15] = -3/(opq1*opq1*opq1*opq1)*x[2]*x[2] - 1/(user->epsilon*user->epsilon);
+
+  ierr = MatSetValues(*B,4,idxm,4,idxm,v,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  if (*J != *B) {
+    ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+    ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  }
+  *flag = SAME_NONZERO_PATTERN;
+  PetscFunctionReturn(0);
+}
+
+/*
+RHS Jacobian (P)
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormJacobianP"
+PetscErrorCode FormJacobianP(TS ts,PetscReal t,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ctx)
+{
+  User              *user = (User*)ctx;
+  PetscErrorCode    ierr;
+  PetscScalar       v[16],*x,opq1;
+  PetscInt          idxm[4] = {0,1,2,3};
+
+  PetscFunctionBeginUser;
+  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
+  opq1 = 1 + x[1];
+ v[0]  = 0;                       v[1]  = 0; v[2] = 0;                           v[3]  = 0;
+ v[4]  = 0;                       v[5]  = 1; v[6] = 0;                           v[5]  = 0;
+ v[8]  = 0;                       v[9]  = 0; v[10] = -opq1*PetscCosScalar(x[0]); v[9]  = -PetscSinScalar(x[0]);
+ v[12] = 2/(opq1*opq1*opq1)*x[2]; v[13] = 0; v[14] = -PetscSinScalar(x[0]);      v[15] = -3/(opq1*opq1*opq1*opq1)*x[2]*x[2] - 1/(user->epsilon*user->epsilon);
+
+  ierr = MatSetValues(*B,4,idxm,4,idxm,v,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
+  ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  if (*J != *B) {
+    ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+    ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  }
+  *flag = SAME_NONZERO_PATTERN;
+  PetscFunctionReturn(0);
+}

# File src/ts/examples/tutorials/makefile

• Ignore whitespace
 LOCDIR          = src/ts/examples/tutorials/
 EXAMPLESC       = ex1.c ex2.c ex3.c ex4.c ex5.c ex6.c ex7.c ex8.c \
                 ex9.c ex10.c ex12.c ex13.c ex14.c ex15.c ex16.c ex17.c \
-                ex19.c ex20.c ex21.c ex22.c ex23.c ex24.c ex25.c ex26.c ex28.c ex30.cxx
+                ex19.c ex20.c ex21.c ex22.c ex23.c ex24.c ex25.c ex26.c \
+                ex28.c ex30.cxx  ex32.c ex33.c
 EXAMPLESF       = ex1f.F ex2f.F ex22f.F # ex22f_mf.F90
 EXAMPLESFH      = ex2f.h
 MANSEC          = TS
 	-${CLINKER} -o ex30 ex30.o${PETSC_TS_LIB}
 	${RM} ex30.o   +ex32: ex32.o chkopts + -${CLINKER} -o ex32 ex32.o ${PETSC_TS_LIB} +${RM} ex32.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/ex30_2.out ex30_2.tmp || printf "${PWD}\nPossible problem with ex30_2, diffs above \n=========================================\n"; \
 	   ${RM} -f ex30_2.tmp   +runex32: + -@${MPIEXEC} -n 1 ./ex32 -ts_monitor_lg_solution > ex32.tmp 2>&1;     \
+    ${DIFF} output/ex32.out ex32.tmp || printf "${PWD}\nPossible problem with ex32, diffs above \n========================================\n"; \
+    ${RM} -f ex32.tmp + +runex33: + -@${MPIEXEC} -n 1 ./ex33 -ts_monitor_lg_solution > ex33.tmp 2>&1;     \
+    ${DIFF} output/ex33.out ex33.tmp || printf "${PWD}\nPossible problem with ex33, diffs above \n========================================\n"; \
+    ${RM} -f ex33.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 \  ex10.PETSc runex10 runex10_2 runex10_3 ex10.rm \  ex12.PETSc ex12.rm ex13.PETSc runex13 runex13_2 runex13_3 ex13.rm\  ex15.PETSc runex15 runex15_2 runex15_3 runex15_4 runex15_5 ex15.rm \ - ex17.PETSc runex17 runex17_2 ex17.rm ex22.PETSc ex22.rm + ex17.PETSc runex17 runex17_2 ex17.rm ex22.PETSc ex22.rm \ + ex32.PETSc runex32 ex32.rm \ + ex33.PETSc runex33 ex33.rm  TESTEXAMPLES_C_NOCOMPLEX = ex9.PETSc runex9 runex9_2 runex9_3 ex9.rm  TESTEXAMPLES_C_X =  TESTEXAMPLES_FORTRAN = ex1f.PETSc runex1f ex1f.rm ex2f.PETSc runex2f ex2f.rm ex22f.PETSc ex22f.rm # ex22f_mf.PETSc ex22f_mf.rm # File src/ts/examples/tutorials/output/ex32.out • Ignore whitespace +CONVERGED_TIME at time 50 after 500 steps # File src/ts/examples/tutorials/output/ex33.out • Ignore whitespace +CONVERGED_TIME at time 10.05 after 67 steps # File src/ts/impls/makefile • Ignore whitespace    ALL: lib   -DIRS = explicit implicit pseudo python arkimex rosw eimex +DIRS = explicit implicit pseudo python arkimex rosw eimex symp  LOCDIR = src/ts/impls/  MANSEC = TS   # File src/ts/impls/symp/makefile • Ignore whitespace + +ALL: lib + +MANSEC = TS +LOCDIR = src/ts/impls/symp/ +DIRS = sympeuler + +include${PETSC_DIR}/conf/variables
+include ${PETSC_DIR}/conf/rules +include${PETSC_DIR}/conf/test
+

# File src/ts/impls/symp/sympeuler/makefile

• Ignore whitespace
+
+ALL: lib
+
+CFLAGS   =
+FFLAGS   =
+SOURCEC  = sympeuler.c
+SOURCEF  =
+SOURCEH  =
+LIBBASE  = libpetscts
+MANSEC   = TS
+LOCDIR   = src/ts/impls/symp/sympeuler/
+
+include ${PETSC_DIR}/conf/variables +include${PETSC_DIR}/conf/rules
+include \${PETSC_DIR}/conf/test
+

# File src/ts/impls/symp/sympeuler/sympeuler.c

• Ignore whitespace
+#include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
+
+typedef struct {
+  Vec update;  /* Work vector for explicit steps (one of the work vectors below could be used instead)*/
+  Vec X0,Xdot; /* work vectors for implicit steps */
+  PetscReal stage_time;
+} TS_SympEuler;
+
+#undef __FUNCT__
+#define __FUNCT__ "TSStep_SympEuler"
+static PetscErrorCode TSStep_SympEuler(TS ts)
+{
+
+  PetscErrorCode      ierr;
+  TS_SympEuler        *sympeuler = (TS_SympEuler*)ts->data;
+  Vec                 sol = ts->vec_sol,update = sympeuler->update;
+  PetscInt            its,lits;
+  SNESConvergedReason snesreason;
+
+  PetscFunctionBegin;
+
+  ierr = TSPreStep(ts);CHKERRQ(ierr);
+  ierr = TSPreStage(ts,ts->ptime);CHKERRQ(ierr);
+
+  sympeuler->stage_time = ts->ptime + ts->time_step; /* full time step (no adaptation) */
+
+  /*  Implicit Step in P */
+  /* Note that for separable systems, this can (and should) just be another euler step. */
+  ierr = VecCopy(sol,sympeuler->X0);CHKERRQ(ierr);
+  ierr = SNESSolve(ts->snes,NULL,sol);CHKERRQ(ierr);
+  ierr = SNESGetIterationNumber(ts->snes,&its);CHKERRQ(ierr);
+  ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr);
+  ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr);
+  ts->snes_its += its; ts->ksp_its += lits;
+
+  /* Explicit Step in Q */
+  /* If the user specified a nontrivial IFunction, this implicit step isn't going to pick it up */
+  ierr = TSComputeRHSPartitionFunction(ts,TS_SYMP_PARTITION,TS_SYMP_Q_SLOT,ts->ptime,sol,update);CHKERRQ(ierr);
+  ierr = VecAXPY(sol,ts->time_step,update);CHKERRQ(ierr);
+  ts->ptime += ts->time_step;
+  ts->steps++;
+
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SNESTSFormFunction_SympEuler"
+static PetscErrorCode SNESTSFormFunction_SympEuler(SNES snes,Vec x,Vec y,TS ts)
+{
+  TS_SympEuler    *sympeuler = (TS_SympEuler*)ts->data;
+  PetscErrorCode  ierr;
+  DM              dm,dmsave;
+  Vec             X0 = sympeuler->X0,Xdot = sympeuler->Xdot;
+  PetscReal       shift = 1./(ts->time_step);
+
+  PetscFunctionBegin;
+  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
+  ierr = VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);CHKERRQ(ierr);
+
+  /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
+  dmsave = ts->dm;
+  ts->dm = dm;
+
+  /* Note that the P step being implicit is hard-coded */
+  ierr   = TSComputeIPartitionFunction(ts,TS_SYMP_PARTITION,TS_SYMP_P_SLOT,sympeuler->stage_time,x,Xdot,y,PETSC_FALSE);CHKERRQ(ierr);
+  ts->dm = dmsave;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SNESTSFormJacobian_SympEuler"
+static PetscErrorCode SNESTSFormJacobian_SympEuler(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *str,TS ts)
+{
+  TS_SympEuler   *sympeuler  = (TS_SympEuler*)ts->data;
+  PetscErrorCode ierr;
+  DM             dm,dmsave;
+  PetscReal      shift = 1./(ts->time_step);
+  Vec            Xdot = sympeuler->Xdot;
+
+  PetscFunctionBegin;
+  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
+
+  /* th->Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
+  dmsave = ts->dm;
+  ts->dm = dm;
+
+  /* Note that the P step being implicit is hard-coded. This can/should change at some point to allow both sympeuler variants */
+  ierr   = TSComputeIPartitionJacobian(ts,TS_SYMP_PARTITION,TS_SYMP_P_SLOT,sympeuler->stage_time,x,Xdot,shift,A,B,str,PETSC_FALSE);CHKERRQ(ierr);
+  ts->dm = dmsave;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "TSSetUp_SympEuler"
+static PetscErrorCode TSSetUp_SympEuler(TS ts)
+{
+  TS_SympEuler       *sympeuler = (TS_SympEuler*)ts->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecDuplicate(ts->vec_sol,&sympeuler->update);CHKERRQ(ierr);
+  ierr = VecDuplicate(ts->vec_sol,&sympeuler->X0);CHKERRQ(ierr);
+  ierr = VecDuplicate(ts->vec_sol,&sympeuler->Xdot);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSReset_SympEuler"
+static PetscErrorCode TSReset_SympEuler(TS ts)
+{
+  TS_SympEuler       *sympeuler = (TS_SympEuler*)ts->data;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecDestroy(&sympeuler->update);CHKERRQ(ierr);
+  ierr = VecDestroy(&sympeuler->X0);CHKERRQ(ierr);
+  ierr = VecDestroy(&sympeuler->Xdot);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSDestroy_SympEuler"
+static PetscErrorCode TSDestroy_SympEuler(TS ts)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = TSReset_SympEuler(ts);CHKERRQ(ierr);
+  ierr = PetscFree(ts->data);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+/*------------------------------------------------------------*/
+
+#undef __FUNCT__
+#define __FUNCT__ "TSSetFromOptions_SympEuler"
+static PetscErrorCode TSSetFromOptions_SympEuler(TS ts)
+{
+  PetscFunctionBegin;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSView_SympEuler"
+static PetscErrorCode TSView_SympEuler(TS ts,PetscViewer viewer)
+{
+  PetscFunctionBegin;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSInterpolate_SympEuler"
+static PetscErrorCode TSInterpolate_SympEuler(TS ts,PetscReal t,Vec X)
+{
+  PetscReal      alpha = (ts->ptime - t)/ts->time_step;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecAXPBY(ts->vec_sol,1.0-alpha,alpha,X);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+/* ------------------------------------------------------------ */
+
+/*MC
+  TSSYMPEULER - ODE solver which takes Backwards Euler steps in the 'P' variables and Euler steps in the 'Q' variables
+
+  Level: beginner
+
+.seealso:  TSCreate(), TS, TSSetType(), TSEULER, TSBEULER
+
+M*/
+#undef __FUNCT__
+#define __FUNCT__ "TSCreate_SympEuler"
+PETSC_EXTERN PetscErrorCode TSCreate_SympEuler(TS ts)
+{
+  TS_SympEuler       *sympeuler;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ts->ops->setup           = TSSetUp_SympEuler;
+  ts->ops->step            = TSStep_SympEuler;
+  ts->ops->reset           = TSReset_SympEuler;
+  ts->ops->destroy         = TSDestroy_SympEuler;
+  ts->ops->setfromoptions  = TSSetFromOptions_SympEuler;
+  ts->ops->view            = TSView_SympEuler;
+  ts->ops->interpolate     = TSInterpolate_SympEuler;
+  ts->ops->snesfunction    = SNESTSFormFunction_SympEuler;
+  ts->ops->snesjacobian    = SNESTSFormJacobian_SympEuler;
+
+  ierr = PetscNewLog(ts,TS_SympEuler,&sympeuler);CHKERRQ(ierr);
+  ts->data = (void*)sympeuler;
+  PetscFunctionReturn(0);
+}

# File src/ts/interface/ts.c

• Ignore whitespace
 .  ts - the TS context obtained from TSCreate()

    Options Database Keys:
-+  -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSGL, TSSSP
++  -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSGL, TSSSP, TSEIMEX, TSARKIMEX, TSSYMPEULER
 .  -ts_max_steps maxsteps - maximum number of time-steps to take
 .  -ts_final_time time - maximum time to compute to
 .  -ts_dt dt - initial time step

# File src/ts/interface/tsregall.c

• Ignore whitespace
 PETSC_EXTERN PetscErrorCode TSCreate_ARKIMEX(TS);
 PETSC_EXTERN PetscErrorCode TSCreate_RosW(TS);
 PETSC_EXTERN PetscErrorCode TSCreate_EIMEX(TS);
+PETSC_EXTERN PetscErrorCode TSCreate_SympEuler(TS);

 #undef __FUNCT__
 #define __FUNCT__ "TSRegisterAll"
   ierr = TSRegister(TSARKIMEX,  TSCreate_ARKIMEX);CHKERRQ(ierr);
   ierr = TSRegister(TSROSW,     TSCreate_RosW);CHKERRQ(ierr);
   ierr = TSRegister(TSEIMEX,    TSCreate_EIMEX);CHKERRQ(ierr);
+  ierr = TSRegister(TSSYMPEULER,TSCreate_SympEuler); CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }