1 #include <petscts.h> 2 #include <stdio.h> 3 4 #define NEW_VERSION // Applicable for the new features; avoid this for the old (current) TSEvent code 5 6 static char help[] = "Simple linear problem with events\n" 7 "x_dot = 0.2*y\n" 8 "y_dot = -0.2*x\n" 9 "Using two event functions = piecewise-polynomials, zeros = 1 (rank-0), 9 (last rank)\n" 10 "Options:\n" 11 "-dir d : zero-crossing direction for events\n" 12 "-flg : additional output in Postevent\n" 13 "-restart : flag for TSRestartStep() in PostEvent\n" 14 "-dtpost x : if x > 0, then on even PostEvent calls dt_postevent = x is set, on odd PostEvent calls dt_postevent = 0 is set,\n" 15 " if x == 0, nothing happens\n"; 16 17 #define MAX_NFUNC 100 // max event functions per rank 18 #define MAX_NEV 5000 // max zero crossings for each rank 19 20 typedef struct { 21 PetscMPIInt rank, size; 22 PetscReal pi; 23 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 24 PetscReal evres[MAX_NEV]; // times of found zero-crossings 25 PetscInt evnum[MAX_NEV]; // number of zero-crossings at each time 26 PetscInt cnt; // counter 27 PetscBool flg; // flag for additional print in PostEvent 28 PetscBool restart; // flag for TSRestartStep() in PostEvent 29 PetscReal dtpost; // post-event step 30 PetscInt postcnt; // counter for PostEvent calls 31 } AppCtx; 32 33 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 34 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 35 36 int main(int argc, char **argv) 37 { 38 TS ts; 39 Mat A; 40 Vec sol; 41 PetscInt n, dir0, m = 0; 42 PetscInt dir[MAX_NFUNC], inds[2]; 43 PetscBool term[MAX_NFUNC]; 44 PetscScalar *x, vals[4]; 45 AppCtx ctx; 46 47 PetscFunctionBeginUser; 48 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 49 setbuf(stdout, NULL); 50 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 51 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 52 ctx.pi = PetscAcosReal(-1.0); 53 ctx.cnt = 0; 54 ctx.flg = PETSC_FALSE; 55 ctx.restart = PETSC_FALSE; 56 ctx.dtpost = 0; 57 ctx.postcnt = 0; 58 59 // The linear problem has a 2*2 matrix. The matrix is constant 60 if (ctx.rank == 0) m = 2; 61 inds[0] = 0; 62 inds[1] = 1; 63 vals[0] = 0; 64 vals[1] = 0.2; 65 vals[2] = -0.2; 66 vals[3] = 0; 67 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A)); 68 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 69 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 70 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 71 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 72 73 PetscCall(MatCreateVecs(A, &sol, NULL)); 74 PetscCall(VecGetArray(sol, &x)); 75 if (ctx.rank == 0) { // initial conditions 76 x[0] = 0; // sin(0) 77 x[1] = 1; // cos(0) 78 } 79 PetscCall(VecRestoreArray(sol, &x)); 80 81 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 82 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 83 84 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 85 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL)); 86 87 PetscCall(TSSetTimeStep(ts, 0.1)); 88 PetscCall(TSSetType(ts, TSBEULER)); 89 PetscCall(TSSetMaxSteps(ts, 10000)); 90 PetscCall(TSSetMaxTime(ts, 10.0)); 91 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 92 PetscCall(TSSetFromOptions(ts)); 93 94 // Set the event handling 95 dir0 = 0; 96 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction 97 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 98 PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep() 99 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step 100 101 n = 0; // event counter 102 if (ctx.rank == 0) { // first event -- on rank-0 103 dir[n] = dir0; 104 term[n++] = PETSC_FALSE; 105 } 106 if (ctx.rank == ctx.size - 1) { // second event -- on last rank 107 dir[n] = dir0; 108 term[n++] = PETSC_FALSE; 109 } 110 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 111 112 // Solution 113 PetscCall(TSSolve(ts, sol)); 114 115 // The 3 columns printed are: [RANK] [num. of events at the given time] [time of event] 116 for (PetscInt j = 0; j < ctx.cnt; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d\t%" PetscInt_FMT "\t%g\n", ctx.rank, ctx.evnum[j], (double)ctx.evres[j])); 117 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 118 119 PetscCall(MatDestroy(&A)); 120 PetscCall(TSDestroy(&ts)); 121 PetscCall(VecDestroy(&sol)); 122 123 PetscCall(PetscFinalize()); 124 return 0; 125 } 126 127 /* 128 User callback for defining the event-functions 129 */ 130 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 131 { 132 PetscInt n = 0; 133 AppCtx *Ctx = (AppCtx *)ctx; 134 135 PetscFunctionBeginUser; 136 // for the test purposes, event-functions are defined based on t 137 // first event -- on rank-0 138 if (Ctx->rank == 0) { 139 if (t < 2.0) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.0, 12)); 140 else gval[n++] = 0.5; 141 } 142 143 // second event -- on last rank 144 if (Ctx->rank == Ctx->size - 1) { 145 if (t > 8.0) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.0, 12)); 146 else gval[n++] = 0.25; 147 } 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150 151 /* 152 User callback for the post-event stuff 153 */ 154 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 155 { 156 AppCtx *Ctx = (AppCtx *)ctx; 157 158 PetscFunctionBeginUser; 159 if (Ctx->flg) { 160 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 161 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 162 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 163 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 164 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 165 } 166 167 if (Ctx->cnt < MAX_NEV && nev_zero > 0) { 168 Ctx->evres[Ctx->cnt] = t; 169 Ctx->evnum[Ctx->cnt++] = nev_zero; 170 } 171 172 #ifdef NEW_VERSION 173 Ctx->postcnt++; // sync 174 if (Ctx->dtpost > 0) { 175 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 176 else PetscCall(TSSetPostEventStep(ts, 0)); 177 } 178 #endif 179 180 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 181 PetscFunctionReturn(PETSC_SUCCESS); 182 } 183 /*---------------------------------------------------------------------------------------------*/ 184 /*TEST 185 test: 186 suffix: 0s1 187 output_file: output/ex2_0s1.out 188 args: -dir 0 189 args: -restart {{0 1}} 190 args: -dtpost {{0 0.25}} 191 args: -ts_event_post_event_step {{0 0.31}} 192 args: -ts_type {{beuler rk}} 193 args: -ts_adapt_type {{none basic}} 194 nsize: 1 195 196 test: 197 suffix: 0s4 198 output_file: output/ex2_0s4.out 199 args: -dir 0 200 args: -restart {{0 1}} 201 args: -dtpost {{0 0.25}} 202 args: -ts_event_post_event_step {{0 0.31}} 203 args: -ts_type {{beuler rk}} 204 args: -ts_adapt_type {{none basic}} 205 nsize: 4 206 filter: sort 207 filter_output: sort 208 209 test: 210 suffix: pos 211 output_file: output/ex2_pos.out 212 args: -dir 1 213 args: -restart {{0 1}} 214 args: -dtpost {{0 0.25}} 215 args: -ts_event_post_event_step {{0 0.31}} 216 args: -ts_type {{beuler rk}} 217 args: -ts_adapt_type {{none basic}} 218 nsize: {{1 4}} 219 filter: sort 220 filter_output: sort 221 222 test: 223 suffix: ns1 224 output_file: output/ex2_ns1.out 225 args: -dir -1 226 args: -restart {{0 1}} 227 args: -dtpost {{0 0.25}} 228 args: -ts_event_post_event_step {{0 0.31}} 229 args: -ts_type {{beuler rk}} 230 args: -ts_adapt_type {{none basic}} 231 nsize: 1 232 233 test: 234 suffix: ns4 235 output_file: output/ex2_ns4.out 236 args: -dir -1 237 args: -restart {{0 1}} 238 args: -dtpost {{0 0.25}} 239 args: -ts_event_post_event_step {{0 0.31}} 240 args: -ts_type {{beuler rk}} 241 args: -ts_adapt_type {{none basic}} 242 nsize: 4 243 filter: sort 244 filter_output: sort 245 TEST*/ 246