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 */ 8 9 #include <petscts.h> 10 11 static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 12 static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *); 13 14 static PetscErrorCode Event(TS, PetscReal, Vec, PetscReal *, void *); 15 static PetscErrorCode PostEvent(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *); 16 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; 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 27 28 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 29 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)); 36 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)); 48 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)); 55 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)); 70 71 PetscCall(TSSolve(ts, NULL)); 72 73 PetscCall(TSDestroy(&ts)); 74 PetscCall(PetscFinalize()); 75 return 0; 76 } 77 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 } 84 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 } 92 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 } 101 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; 106 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 } 114 115 /*TEST 116 117 test: 118 args: -ts_type beuler -ts_dt 0.1 -ts_event_monitor 119 120 test: 121 suffix: 2 122 args: -ts_type beuler -ts_dt 0.2 -ts_event_monitor 123 124 test: 125 suffix: 3 126 args: -ts_type beuler -ts_dt 0.5 -ts_event_monitor 127 TEST*/ 128