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

  9: #include <petscts.h>

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

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

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

 25:   PetscFunctionBeginUser;
 26:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 28:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));

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

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

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

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

 71:   PetscCall(TSSolve(ts, NULL));

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

 78: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
 79: {
 80:   PetscFunctionBeginUser;
 81:   PetscCall(VecSet(f, (PetscReal)1));
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
 86: {
 87:   PetscFunctionBeginUser;
 88:   PetscCall(MatZeroEntries(B));
 89:   if (B != A) PetscCall(MatZeroEntries(A));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

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

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

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

115: /*TEST

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

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

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