1 #include <petscts.h> 2 #include <stdio.h> 3 4 #define NEW_VERSION // Applicable for the new features; avoid this for the older PETSc versions (without TSSetPostEventStep()) 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 10 "The following event functions are involved:\n" 11 "- two polynomial event functions on rank-0 and last-rank (with zeros: 1.05, 9.05[terminating])\n" 12 "- one event function on rank = '1%size', equal to V*sin(pi*t), zeros = 1,...,10\n" 13 "After each event location the tolerance for the sin() event is multiplied by 4\n" 14 15 "Options:\n" 16 "-dir d : zero-crossing direction for events: 0, 1, -1\n" 17 "-flg : additional output in Postevent\n" 18 "-errtol e : error tolerance, for printing 'pass/fail' for located events (1e-5 by default)\n" 19 "-restart : flag for TSRestartStep() in PostEvent\n" 20 "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n" 21 " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n" 22 " if x == 0, nothing happens\n" 23 "-V {float}: scaling of the sin() event function; for small V this event is triggered by the function values,\n" 24 " for large V the event is triggered by the small step size\n" 25 "-change5 : flag to change the state vector at t=5 PostEvent\n"; 26 27 #define MAX_NFUNC 100 // max event functions per rank 28 #define MAX_NEV 5000 // max zero crossings for each rank 29 30 typedef struct { 31 PetscMPIInt rank, size; 32 PetscReal pi; 33 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 34 PetscReal evres[MAX_NEV]; // times of found zero-crossings 35 PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking 36 PetscInt cnt; // counter 37 PetscInt cntref; // actual length of 'ref' on the given rank 38 PetscBool flg; // flag for additional print in PostEvent 39 PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default) 40 PetscBool restart; // flag for TSRestartStep() in PostEvent 41 PetscReal dtpost; // post-event step 42 PetscInt postcnt; // counter for PostEvent calls 43 PetscReal V; // vertical scaling for sin() 44 PetscReal vtol[MAX_NFUNC]; // vtol array, with extra storage 45 PetscBool change5; // flag to change the state vector at t=5 PostEvent 46 } AppCtx; 47 48 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 49 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 50 51 int main(int argc, char **argv) 52 { 53 TS ts; 54 Mat A; 55 Vec sol; 56 PetscInt n, dir0, m = 0; 57 PetscReal tol = 1e-7; 58 PetscInt dir[MAX_NFUNC], inds[2]; 59 PetscBool term[MAX_NFUNC]; 60 PetscScalar *x, vals[4]; 61 AppCtx ctx; 62 TSConvergedReason reason; 63 64 PetscFunctionBeginUser; 65 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 66 setbuf(stdout, NULL); 67 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 68 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 69 ctx.pi = PetscAcosReal(-1.0); 70 ctx.cnt = 0; 71 ctx.cntref = 0; 72 ctx.flg = PETSC_FALSE; 73 ctx.errtol = 1e-5; 74 ctx.restart = PETSC_FALSE; 75 ctx.dtpost = 0; 76 ctx.postcnt = 0; 77 ctx.V = 1.0; 78 ctx.change5 = PETSC_FALSE; 79 80 // The linear problem has a 2*2 matrix. The matrix is constant 81 if (ctx.rank == 0) m = 2; 82 inds[0] = 0; 83 inds[1] = 1; 84 vals[0] = 0; 85 vals[1] = 0.2; 86 vals[2] = -0.2; 87 vals[3] = 0; 88 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A)); 89 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 90 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 91 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 92 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 93 94 PetscCall(MatCreateVecs(A, &sol, NULL)); 95 PetscCall(VecGetArray(sol, &x)); 96 if (ctx.rank == 0) { // initial conditions 97 x[0] = 0; // sin(0) 98 x[1] = 1; // cos(0) 99 } 100 PetscCall(VecRestoreArray(sol, &x)); 101 102 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 103 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 104 105 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 106 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL)); 107 108 PetscCall(TSSetTimeStep(ts, 0.1)); 109 PetscCall(TSSetType(ts, TSBEULER)); 110 PetscCall(TSSetMaxSteps(ts, 10000)); 111 PetscCall(TSSetMaxTime(ts, 10.0)); 112 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 113 PetscCall(TSSetFromOptions(ts)); 114 115 // Set the event handling 116 dir0 = 0; 117 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction 118 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 119 PetscCall(PetscOptionsGetReal(NULL, NULL, "-errtol", &ctx.errtol, NULL)); // error tolerance for located events 120 PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep() 121 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step 122 PetscCall(PetscOptionsGetReal(NULL, NULL, "-V", &ctx.V, NULL)); 123 PetscCall(PetscOptionsGetBool(NULL, NULL, "-change5", &ctx.change5, NULL)); // flag to change the state vector at t=5 PostEvent 124 125 n = 0; // event counter 126 if (ctx.rank == 0) { // first event -- on rank-0 127 ctx.vtol[n] = tol * 10; 128 dir[n] = dir0; 129 term[n++] = PETSC_FALSE; 130 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.05; 131 } 132 if (ctx.rank == ctx.size - 1) { // second event (with termination) -- on last rank 133 ctx.vtol[n] = tol * 10; 134 dir[n] = dir0; 135 term[n++] = PETSC_TRUE; 136 if (dir0 <= 0) ctx.ref[ctx.cntref++] = 9.05; 137 } 138 if (ctx.rank == 1 % ctx.size) { // third event -- on rank = 1%ctx.size 139 ctx.vtol[n] = tol; 140 dir[n] = dir0; 141 term[n++] = PETSC_FALSE; 142 143 for (PetscInt i = 1; i < MAX_NEV - 2; i++) { 144 if (i % 2 == 1 && dir0 <= 0) ctx.ref[ctx.cntref++] = i; 145 if (i % 2 == 0 && dir0 >= 0) ctx.ref[ctx.cntref++] = i; 146 } 147 } 148 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref)); 149 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 150 PetscCall(TSSetEventTolerances(ts, tol, ctx.vtol)); 151 152 // Solution 153 PetscCall(TSSolve(ts, sol)); 154 PetscCall(TSGetConvergedReason(ts, &reason)); 155 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CONVERGED REASON: %" PetscInt_FMT " (TS_CONVERGED_EVENT == %" PetscInt_FMT ")\n", (PetscInt)reason, (PetscInt)TS_CONVERGED_EVENT)); 156 157 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"] 158 for (PetscInt j = 0; j < ctx.cnt; j++) { 159 PetscReal err = 10.0; 160 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]); 161 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d\t%g\t%g\t%s\n", ctx.rank, (double)ctx.evres[j], (double)err, err < ctx.errtol ? "pass" : "fail")); 162 } 163 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 164 165 PetscCall(MatDestroy(&A)); 166 PetscCall(TSDestroy(&ts)); 167 PetscCall(VecDestroy(&sol)); 168 169 PetscCall(PetscFinalize()); 170 return 0; 171 } 172 173 /* 174 User callback for defining the event-functions 175 */ 176 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 177 { 178 PetscInt n = 0; 179 AppCtx *Ctx = (AppCtx *)ctx; 180 181 PetscFunctionBeginUser; 182 // for the test purposes, event-functions are defined based on t 183 // first event -- on rank-0 184 if (Ctx->rank == 0) { 185 if (t < 2.05) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.05, 12)); 186 else gval[n++] = 0.5; 187 } 188 189 // second event -- on last rank 190 if (Ctx->rank == Ctx->size - 1) { 191 if (t > 8.05) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.05, 12)); 192 else gval[n++] = 0.25; 193 } 194 195 // third event -- on rank = 1%ctx.size 196 if (Ctx->rank == 1 % Ctx->size) { gval[n++] = Ctx->V * PetscSinReal(Ctx->pi * t); } 197 PetscFunctionReturn(PETSC_SUCCESS); 198 } 199 200 /* 201 User callback for the post-event stuff 202 */ 203 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 204 { 205 PetscInt n = 0; 206 PetscScalar *x; 207 AppCtx *Ctx = (AppCtx *)ctx; 208 209 PetscFunctionBeginUser; 210 if (Ctx->flg) { 211 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 212 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 213 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 214 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 215 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 216 } 217 218 // change the state vector near t=5.0 219 if (PetscAbsReal(t - 5.0) < 0.01 && Ctx->change5) { 220 PetscCall(VecGetArray(U, &x)); 221 if (Ctx->rank == 0) x[1] = -x[1]; 222 PetscCall(VecRestoreArray(U, &x)); 223 } 224 225 // update vtol's 226 if (Ctx->rank == 0) n++; // first event -- on rank-0 227 if (Ctx->rank == Ctx->size - 1) n++; // second event -- on last rank 228 if (Ctx->rank == 1 % Ctx->size) { // third event -- on rank = 1%ctx.size 229 if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "vtol for sin: %g -> ", (double)Ctx->vtol[n])); 230 Ctx->vtol[n] *= 4; 231 if (PetscAbsReal(t - 5.0) < 0.01) Ctx->vtol[n] /= 100; // one-off decrease 232 if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g\n", (double)Ctx->vtol[n])); 233 n++; 234 } 235 PetscCall(TSSetEventTolerances(ts, 0, Ctx->vtol)); 236 237 if (Ctx->cnt + nev_zero < MAX_NEV) 238 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing 239 240 #ifdef NEW_VERSION 241 Ctx->postcnt++; // sync 242 if (Ctx->dtpost > 0) { 243 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 244 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE)); 245 } 246 #endif 247 248 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 249 PetscFunctionReturn(PETSC_SUCCESS); 250 } 251 /*---------------------------------------------------------------------------------------------*/ 252 /* 253 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1, 254 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to 255 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed, 256 simply remove this option altogether. This will result in using the defaults 257 (which is PETSC_DECIDE). 258 */ 259 /*TEST 260 test: 261 suffix: V 262 output_file: output/ex3_V.out 263 args: -ts_type beuler 264 args: -ts_adapt_type basic 265 args: -V {{1e2 1e5 1e8}} 266 args: -ts_adapt_dt_min 1e-6 267 args: -change5 {{0 1}} 268 nsize: 1 269 270 test: 271 suffix: neu1 272 output_file: output/ex3_neu1.out 273 args: -dir 0 274 args: -V 1e5 275 args: -ts_adapt_dt_min 1e-6 276 args: -restart 1 277 args: -dtpost 0.24 278 args: -ts_event_post_event_step 0.31 279 args: -ts_type {{beuler rk}} 280 args: -ts_adapt_type {{none basic}} 281 nsize: 1 282 283 test: 284 suffix: neu2 285 output_file: output/ex3_neu2.out 286 args: -dir 0 287 args: -V 1e5 288 args: -ts_adapt_dt_min 1e-6 289 args: -restart 1 290 args: -dtpost 0 291 args: -ts_event_post_event_step {{-1 0.31}} 292 args: -ts_type rk 293 args: -ts_adapt_type {{none basic}} 294 nsize: 2 295 filter: sort 296 filter_output: sort 297 298 test: 299 suffix: neu4 300 output_file: output/ex3_neu4.out 301 args: -dir 0 302 args: -V 1e5 303 args: -ts_adapt_dt_min 1e-6 304 args: -restart {{0 1}} 305 args: -dtpost 0.24 306 args: -ts_event_post_event_step 0.21 307 args: -ts_type beuler 308 args: -ts_adapt_type {{none basic}} 309 nsize: 4 310 filter: sort 311 filter_output: sort 312 313 test: 314 suffix: pos1 315 output_file: output/ex3_pos1.out 316 args: -dir 1 317 args: -V 1e5 318 args: -ts_adapt_dt_min 1e-6 319 args: -restart 0 320 args: -dtpost 0.24 321 args: -ts_type {{beuler rk}} 322 args: -ts_adapt_type {{none basic}} 323 nsize: 1 324 325 test: 326 suffix: pos2 327 output_file: output/ex3_pos2.out 328 args: -dir 1 329 args: -V 1e5 330 args: -ts_adapt_dt_min 1e-6 331 args: -restart 1 332 args: -dtpost {{0 0.24}} 333 args: -ts_type rk 334 args: -ts_adapt_type {{none basic}} 335 nsize: 2 336 filter: sort 337 filter_output: sort 338 339 test: 340 suffix: pos4 341 output_file: output/ex3_pos4.out 342 args: -dir 1 343 args: -V 1e9 344 args: -ts_adapt_dt_min 1e-6 345 args: -restart 0 346 args: -dtpost 0 347 args: -ts_event_post_event_step {{-1 0.32}} 348 args: -ts_type beuler 349 args: -ts_adapt_type {{none basic}} 350 args: -change5 1 351 nsize: 4 352 filter: sort 353 filter_output: sort 354 355 test: 356 suffix: neg1 357 output_file: output/ex3_neg1.out 358 args: -dir -1 359 args: -V 1e5 360 args: -ts_adapt_dt_min 1e-6 361 args: -restart 1 362 args: -dtpost {{0 0.24}} 363 args: -ts_type {{beuler rk}} 364 args: -ts_adapt_type basic 365 nsize: 1 366 367 test: 368 suffix: neg2 369 output_file: output/ex3_neg2.out 370 args: -dir -1 371 args: -V 1e5 372 args: -ts_adapt_dt_min 1e-6 373 args: -restart 0 374 args: -dtpost {{0 0.24}} 375 args: -ts_type rk 376 args: -ts_adapt_type {{none basic}} 377 nsize: 2 378 filter: sort 379 filter_output: sort 380 381 test: 382 suffix: neg4 383 output_file: output/ex3_neg4.out 384 args: -dir -1 385 args: -V 1e5 386 args: -ts_adapt_dt_min 1e-6 387 args: -restart 0 388 args: -dtpost {{0 0.24}} 389 args: -ts_event_post_event_step 0.3 390 args: -ts_type beuler 391 args: -ts_adapt_type {{none basic}} 392 nsize: 4 393 filter: sort 394 filter_output: sort 395 TEST*/ 396