# File include/petscts.h Modified

• Ignore whitespace
• Hide word diff
 #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

 PETSC_EXTERN PetscErrorCode TSGetIJacobian(TS,Mat*,Mat*,TSIJacobian*,void**);

 PETSC_EXTERN PetscErrorCode TSComputeRHSFunctionLinear(TS,PetscReal,Vec,Vec,void*);
+PETSC_EXTERN PetscErrorCode TSComputeRHSPartitionFunctionLinear(TS,TSPartitionType,TSPartitionSlotType,PetscReal,Vec,Vec,void*);
 PETSC_EXTERN PetscErrorCode TSComputeRHSJacobianConstant(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
 PETSC_EXTERN PetscErrorCode TSComputeIFunctionLinear(TS,PetscReal,Vec,Vec,Vec,void*);
 PETSC_EXTERN PetscErrorCode TSComputeIJacobianConstant(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);

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

• Ignore whitespace
• Hide word diff

    q''(t) + q(t) = \epsilon f(q)

-   where \epsilon << 1
+   where \epsilon may or may not be small

    As a first order system, we  have

                    -- P ---   ----------- Q -------------------
                    ------- Fast --------   ------ Slow --------

-  Different potential integration schemes would split the same RHS differently.
+  Different potential integration schemes would split the same RHS differently,
+  perhaps based on the value of \epsilon
+
+  Accepts an extra option -epsilon
 */

 #include <petscts.h>
   TS ts;
   DM dm;
   User user;
-  PetscReal dt = 0.1, maxtime = 1.0, ftime;
+  PetscScalar epsilon = 1;
+  PetscReal dt = 0.1, maxtime = 50, ftime;
   PetscInt maxsteps = 1000, 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 = 0.01;
+  user.epsilon = epsilon;

   /* Create TS */
   ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
   ierr = DMTSSetRHSPartitionFunction(dm,EXPONENTIAL,EXPONENTIAL_FAST,FormRHSFunctionFast,&user);CHKERRQ(ierr);
   ierr = DMTSSetRHSPartitionFunction(dm,EXPONENTIAL,EXPONENTIAL_SLOW,FormRHSFunctionSlow,&user);CHKERRQ(ierr);

-  // Jacobian Functions ...
   // TSSetRHSJacobian does error checking and some matrix wrangling
   //  it also sets something in snes, which may be essential for implicit methods..
-  //ierr = DMTSSetRHSJacobian(dm,FormJacobian,&user);CHKERRQ(ierr);

   ierr = TSSetRHSJacobian(ts, J, J, FormJacobian, &user);CHKERRQ(ierr);
   ierr = DMTSSetRHSPartitionJacobian(dm, SYMPLECTIC,SYMPLECTIC_Q, FormJacobianQ, &user);CHKERRQ(ierr);

   // Note that this still seems shaky - it relies on all the FormJacobian functions being very interchangeable, since the setup that happens in TSSetRHSJacobian only happens once. Also, not sure if that registers the wrong function with the SNES object!

+  /* Set TS Type  */
+  ierr = TSSetType(ts,TSSYMPEULER);CHKERRQ(ierr);
+
   /* Set from options */
   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);


# File src/ts/examples/tutorials/makefile Modified

• Ignore whitespace
• Hide word diff
     ${RM} -f ex31.tmp    runex32: - -@${MPIEXEC} -n 1 ./ex32 -ts_monitor_lg_solution -ts_final_time 30> ex32.tmp 2>&1;     \
+	-@${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


# File src/ts/examples/tutorials/output/ex32.out Modified

• Ignore whitespace
• Hide word diff
-This output file (output/ex32.out) does not contain meaningful information, hence this example is incomplete.
+CONVERGED_TIME at time 50 after 500 steps

# File src/ts/impls/makefile Modified

• Ignore whitespace
• Hide word diff

 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


• Ignore whitespace
• Hide word diff
+
+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 Added • Ignore whitespace • Hide word diff + +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
+

• Ignore whitespace
• Hide word diff
+/*
+  Code for timestepping with the Symplectic Euler method
+*/
+#include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
+
+typedef struct {
+  Vec update; /* Work vector */
+} TS_SympEuler;
+
+#undef __FUNCT__
+#define __FUNCT__ "TSStep_SympEuler"
+static PetscErrorCode TSStep_SympEuler(TS ts)
+{
+
+  TS_SympEuler    *sympeuler = (TS_SympEuler*)ts->data;
+  Vec            sol = ts->vec_sol,update = sympeuler->update;
+  PetscErrorCode ierr;
+  // ..
+  PetscFunctionBegin;
+
+  // Note - there are two variants on symplectic euler ( which of p or q is treated as implicit)
+  //   Here we choose one arbitrarily, but later when we construct a proper subpackage we can
+  //   have two different integrators (of course, practically speaking, the user could just switch which variables are called P and which are called Q)
+
+  // Implicit Step in Q
+  //..
+  // Explicit Step in P
+  //..
+  // DEBUG - take both explicit steps for now
+  // (This works for the mechanical Lagrangian type problems - the trick is to reproduce
+  //    the speed and simplicity for that system when we can, and deal with the actual implicit solve when we need to)
+  //
+  ierr = TSPreStep(ts);CHKERRQ(ierr);
+  ierr = TSPreStage(ts,ts->ptime);CHKERRQ(ierr);
+  ierr = TSComputeRHSPartitionFunction(ts,SYMPLECTIC,SYMPLECTIC_P,ts->ptime,sol,update);CHKERRQ(ierr);
+  ierr = VecAXPY(sol,ts->time_step,update);CHKERRQ(ierr);
+  ierr = TSComputeRHSPartitionFunction(ts,SYMPLECTIC,SYMPLECTIC_Q,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_5
+#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);
+  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);
+  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);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSComputeLinearStability_SympEuler"
+PetscErrorCode TSComputeLinearStability_SympEuler(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
+{
+  //!! This is just copied from euler !!
+  PetscFunctionBegin;
+  *yr = 1.0 + xr;
+  *yi = xi;
+  PetscFunctionReturn(0);
+}
+/* ------------------------------------------------------------ */
+
+/*MC
+      TSSYMPEULER - ODE solver using the explicit forward Euler method
+
+  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->linearstability = TSComputeLinearStability_SympEuler;
+
+  ierr = PetscNewLog(ts,TS_SympEuler,&sympeuler);CHKERRQ(ierr);
+  ts->data = (void*)sympeuler;
+  PetscFunctionReturn(0);
+}

# File src/ts/interface/ts.c Modified

• Ignore whitespace
• Hide word diff
 .  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
 PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y)
 {
   PetscErrorCode ierr;
+
+  //!! Todo is change to just using the partition function, since there's no point in having this extra stack frame
+  PetscFunctionBegin;
+  ierr = TSComputeRHSPartitionFunction(ts,NONE,DEFAULT,t,U,y);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TSComputeRHSPartitionFunction"
+/*@
+   TSComputeRHSPartitionFunction - Evaluates the right-hand-side function for a given partition and 'slot'
+
+   Collective on TS and Vec
+
+   Input Parameters:
++  ts - the TS context
+.  type - the partition type
+.  slot - the partition slot type
+.  t - current time
+-  U - state vector
+
+   Output Parameter:
+.  y - right hand side
+
+   Note:
+   Most users should not need to explicitly call this routine, as it
+   is used internally within the nonlinear solvers.
+
+   Level: developer
+
+.keywords: TS, compute
+
+.seealso: TSSetRHSFunction(), TSComputeIFunction(), TSComputeRHSFunction()
+@*/
+PetscErrorCode TSComputeRHSPartitionFunction(TS ts,TSPartitionType type, TSPartitionSlotType slot, PetscReal t,Vec U,Vec y)
+{
+  PetscErrorCode ierr;
   TSRHSFunction  rhsfunction;
   TSIFunction    ifunction;
   void           *ctx;

   PetscFunctionBegin;
   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
-  PetscValidHeaderSpecific(U,VEC_CLASSID,3);
-  PetscValidHeaderSpecific(y,VEC_CLASSID,4);
+  PetscValidHeaderSpecific(U,VEC_CLASSID,5);
+  PetscValidHeaderSpecific(y,VEC_CLASSID,6);
   ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
-  ierr = DMTSGetRHSFunction(dm,&rhsfunction,&ctx);CHKERRQ(ierr);
+  ierr = DMTSGetRHSPartitionFunction(dm,type,slot,&rhsfunction,&ctx);CHKERRQ(ierr);
   ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr);

   if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()");
   ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
-
 #undef __FUNCT__
 #define __FUNCT__ "TSComputeSolutionFunction"
 /*@

# File src/ts/interface/tsregall.c Modified

• Ignore whitespace
• Hide word diff
 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);
 }