Actual source code: ex16.c


  2: static char help[] = "Solves the trivial ODE du/dt = 1, u(0) = 0. \n\n";
  3: /*
  4:   This example tests TSEvent's capability to handle complicated cases.
  5:   Test 1: an event close to endpoint of a time step should not be detected twice.
  6:   Test 2: two events in the same time step should be detected in the correct order.
  7:   Test 3: three events in the same time step should be detected in the correct order.
  8: */

 10: #include <petscts.h>

 12: static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
 13: static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);

 15: static PetscErrorCode Event(TS, PetscReal, Vec, PetscScalar *, void *);
 16: static PetscErrorCode PostEvent(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *);

 18: int main(int argc, char **argv)
 19: {
 20:   TS              ts;
 21:   const PetscReal t_end = 2.0;
 22:   Vec             x;
 23:   Vec             f;
 24:   Mat             A;

 27:   PetscInitialize(&argc, &argv, (char *)0, help);

 29:   TSCreate(PETSC_COMM_WORLD, &ts);

 31:   VecCreate(PETSC_COMM_WORLD, &f);
 32:   VecSetSizes(f, 1, PETSC_DECIDE);
 33:   VecSetFromOptions(f);
 34:   VecSetUp(f);
 35:   TSSetRHSFunction(ts, f, RHSFunction, NULL);
 36:   VecDestroy(&f);

 38:   MatCreate(PETSC_COMM_WORLD, &A);
 39:   MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE);
 40:   MatSetFromOptions(A);
 41:   MatSetUp(A);
 42:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 43:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 44:   /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
 45:   MatShift(A, (PetscReal)1);
 46:   MatShift(A, (PetscReal)-1);
 47:   TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL);
 48:   MatDestroy(&A);

 50:   VecCreate(PETSC_COMM_WORLD, &x);
 51:   VecSetSizes(x, 1, PETSC_DECIDE);
 52:   VecSetFromOptions(x);
 53:   VecSetUp(x);
 54:   TSSetSolution(ts, x);
 55:   VecDestroy(&x);

 57:   {
 58:     PetscInt  direction[3];
 59:     PetscBool terminate[3];
 60:     direction[0] = +1;
 61:     terminate[0] = PETSC_FALSE;
 62:     direction[1] = -1;
 63:     terminate[1] = PETSC_FALSE;
 64:     direction[2] = 0;
 65:     terminate[2] = PETSC_FALSE;
 66:     TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL);
 67:   }
 68:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
 69:   TSSetMaxTime(ts, t_end);
 70:   TSSetFromOptions(ts);

 72:   TSSolve(ts, NULL);

 74:   TSDestroy(&ts);
 75:   PetscFinalize();
 76:   return 0;
 77: }

 79: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
 80: {
 82:   VecSet(f, (PetscReal)1);
 83:   return 0;
 84: }

 86: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
 87: {
 89:   MatZeroEntries(B);
 90:   if (B != A) MatZeroEntries(A);
 91:   return 0;
 92: }

 94: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
 95: {
 97:   fvalue[0] = t - 1.1;
 98:   fvalue[1] = 1.2 - t;
 99:   fvalue[2] = t - 1.3;
100:   return 0;
101: }

103: PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx)
104: {
105:   PetscInt           i;
106:   const PetscScalar *a;

109:   TSGetStepNumber(ts, &i);
110:   VecGetArrayRead(x, &a);
111:   PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0]));
112:   VecRestoreArrayRead(x, &a);
113:   return 0;
114: }

116: /*TEST

118:     test:
119:       args: -ts_type beuler -ts_dt 0.1 -ts_event_monitor

121:     test:
122:       suffix: 2
123:       args: -ts_type beuler -ts_dt 0.2 -ts_event_monitor

125:     test:
126:       suffix: 3
127:       args: -ts_type beuler -ts_dt 0.5 -ts_event_monitor
128: TEST*/