1. petsc
  2. PETSc
  3. petsc

Commits

Patrick Sanan  committed ee68e94

new sympletic euler TS, fully explicit

  • Participants
  • Parent commits 4f44656
  • Branches psanan/ts-tests

Comments (0)

Files changed (10)

File include/petscts.h

View file
 #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

View file
 
    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

View file
     ${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

View file
-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

View file
 
 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

View file
+
+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

View file
+
+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

View file
+/*
+  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

View file
 .  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

View file
 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);
 }