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;

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

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

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

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

 50:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 51:   PetscCall(VecSetSizes(x, 1, PETSC_DECIDE));
 52:   PetscCall(VecSetFromOptions(x));
 53:   PetscCall(VecSetUp(x));
 54:   PetscCall(TSSetSolution(ts, x));
 55:   PetscCall(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:     PetscCall(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL));
 67:   }
 68:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
 69:   PetscCall(TSSetMaxTime(ts, t_end));
 70:   PetscCall(TSSetFromOptions(ts));

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

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

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

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

 94: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
 95: {
 96:   PetscFunctionBeginUser;
 97:   fvalue[0] = t - 1.1;
 98:   fvalue[1] = 1.2 - t;
 99:   fvalue[2] = t - 1.3;
100:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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