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 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 sin(pi*t), zeros = 1,...,10\n" 13 "TimeSpan = [0.01, 0.21, 1.01, ..., 6.21, 6.99, 7.21,... 9.21] plus the points: {3, 4, 4+D, 5-D, 5, 6-D, 6, 6+D} with user-defined D\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 "-restart : flag for TSRestartStep() in PostEvent\n" 19 "-term : flag to terminate at 9.05 event (true by default)\n" 20 "-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" 21 " if x == 0, nothing happens\n" 22 "-D z : a small number to define additional TimeSpan points\n"; 23 24 #define MAX_NFUNC 100 // max event functions per rank 25 #define MAX_NEV 5000 // max zero crossings for each rank 26 27 typedef struct { 28 PetscMPIInt rank, size; 29 PetscReal pi; 30 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 31 PetscReal evres[MAX_NEV]; // times of found zero-crossings 32 PetscInt evnum[MAX_NEV]; // number of zero-crossings at each time 33 PetscInt cnt; // counter 34 PetscBool flg; // flag for additional print in PostEvent 35 PetscBool restart; // flag for TSRestartStep() in PostEvent 36 PetscBool term; // flag to terminate at 9.05 event 37 PetscReal dtpost; // post-event step 38 PetscInt postcnt; // counter for PostEvent calls 39 } AppCtx; 40 41 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 42 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 43 44 int main(int argc, char **argv) 45 { 46 TS ts; 47 Mat A; 48 Vec sol; 49 PetscInt n, dir0, m = 0; 50 PetscReal tol = 1e-7, D = 0.02; 51 PetscInt dir[MAX_NFUNC], inds[2]; 52 PetscBool term[MAX_NFUNC], match; 53 PetscScalar *x, vals[4]; 54 PetscReal tspan[28], dtlast; 55 AppCtx ctx; 56 TSConvergedReason reason; 57 TSAdapt adapt; 58 59 PetscFunctionBeginUser; 60 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 61 setbuf(stdout, NULL); 62 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 63 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 64 ctx.pi = PetscAcosReal(-1.0); 65 ctx.cnt = 0; 66 ctx.flg = PETSC_FALSE; 67 ctx.restart = PETSC_FALSE; 68 ctx.term = PETSC_TRUE; 69 ctx.dtpost = 0; 70 ctx.postcnt = 0; 71 72 // The linear problem has a 2*2 matrix. The matrix is constant 73 if (ctx.rank == 0) m = 2; 74 inds[0] = 0; 75 inds[1] = 1; 76 vals[0] = 0; 77 vals[1] = 0.2; 78 vals[2] = -0.2; 79 vals[3] = 0; 80 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A)); 81 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 82 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 83 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 84 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 85 86 PetscCall(MatCreateVecs(A, &sol, NULL)); 87 PetscCall(VecGetArray(sol, &x)); 88 if (ctx.rank == 0) { // initial conditions 89 x[0] = 0; // sin(0) 90 x[1] = 1; // cos(0) 91 } 92 PetscCall(VecRestoreArray(sol, &x)); 93 94 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 95 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 96 97 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 98 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL)); 99 100 PetscCall(TSSetTimeStep(ts, 0.099)); 101 PetscCall(TSSetType(ts, TSBEULER)); 102 PetscCall(TSSetMaxSteps(ts, 10000)); 103 PetscCall(TSSetMaxTime(ts, 10.0)); 104 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 105 106 // Set the event handling 107 dir0 = 0; 108 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction 109 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 110 PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep() 111 PetscCall(PetscOptionsGetBool(NULL, NULL, "-term", &ctx.term, NULL)); // flag to terminate at 9.05 event 112 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step 113 PetscCall(PetscOptionsGetReal(NULL, NULL, "-D", &D, NULL)); // small number for tspan 114 115 n = 0; // event counter 116 if (ctx.rank == 0) { // first event -- on rank-0 117 dir[n] = dir0; 118 term[n++] = PETSC_FALSE; 119 } 120 if (ctx.rank == ctx.size - 1) { // second event (with optional termination) -- on last rank 121 dir[n] = dir0; 122 term[n++] = ctx.term; 123 } 124 if (ctx.rank == 1 % ctx.size) { // third event -- on rank = 1%ctx.size 125 dir[n] = dir0; 126 term[n++] = PETSC_FALSE; 127 } 128 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 129 PetscCall(TSSetEventTolerances(ts, tol, NULL)); 130 131 // Set the time span 132 for (PetscInt i = 0; i < 10; i++) { 133 tspan[2 * i] = 0.01 + i + (i == 7 ? -0.02 : 0); 134 tspan[2 * i + 1] = 0.21 + i; 135 } 136 tspan[20] = 3; 137 tspan[21] = 4; 138 tspan[22] = 4 + D; 139 tspan[23] = 5 - D; 140 tspan[24] = 5; 141 tspan[25] = 6 - D; 142 tspan[26] = 6; 143 tspan[27] = 6 + D; 144 PetscCall(PetscSortReal(28, tspan)); 145 PetscCall(TSSetTimeSpan(ts, 28, tspan)); 146 PetscCall(TSSetFromOptions(ts)); 147 148 // Solution 149 PetscCall(TSSolve(ts, sol)); 150 PetscCall(TSGetConvergedReason(ts, &reason)); 151 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CONVERGED REASON: %" PetscInt_FMT " (TS_CONVERGED_EVENT == %" PetscInt_FMT ")\n", (PetscInt)reason, (PetscInt)TS_CONVERGED_EVENT)); 152 153 // The 3 columns printed are: [RANK] [num. of events at the given time] [time of event] 154 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])); 155 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 156 157 // print the final time step 158 PetscCall(TSGetTimeStep(ts, &dtlast)); 159 PetscCall(TSGetAdapt(ts, &adapt)); 160 PetscCall(PetscObjectTypeCompare((PetscObject)adapt, TSADAPTNONE, &match)); 161 if (match) { 162 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Adapt = none\n")); 163 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Last dt = %g\n", (double)dtlast)); 164 } 165 166 PetscCall(MatDestroy(&A)); 167 PetscCall(TSDestroy(&ts)); 168 PetscCall(VecDestroy(&sol)); 169 170 PetscCall(PetscFinalize()); 171 return 0; 172 } 173 174 /* 175 User callback for defining the event-functions 176 */ 177 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 178 { 179 PetscInt n = 0; 180 AppCtx *Ctx = (AppCtx *)ctx; 181 182 PetscFunctionBeginUser; 183 // for the test purposes, event-functions are defined based on t 184 // first event -- on rank-0 185 if (Ctx->rank == 0) { 186 if (t < 2.05) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.05, 12)); 187 else gval[n++] = 0.5; 188 } 189 190 // second event -- on last rank 191 if (Ctx->rank == Ctx->size - 1) { 192 if (t > 8.05) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.05, 12)); 193 else gval[n++] = 0.25; 194 } 195 196 // third event -- on rank = 1%ctx.size 197 if (Ctx->rank == 1 % Ctx->size) { gval[n++] = PetscSinReal(Ctx->pi * t); } 198 PetscFunctionReturn(PETSC_SUCCESS); 199 } 200 201 /* 202 User callback for the post-event stuff 203 */ 204 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 205 { 206 AppCtx *Ctx = (AppCtx *)ctx; 207 208 PetscFunctionBeginUser; 209 if (Ctx->flg) { 210 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 211 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 212 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 213 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 214 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 215 } 216 217 if (Ctx->cnt < MAX_NEV && nev_zero > 0) { 218 Ctx->evres[Ctx->cnt] = t; 219 Ctx->evnum[Ctx->cnt++] = nev_zero; 220 } 221 222 #ifdef NEW_VERSION 223 Ctx->postcnt++; // sync 224 if (Ctx->dtpost > 0) { 225 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 226 else PetscCall(TSSetPostEventStep(ts, 0)); 227 } 228 #endif 229 230 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 231 PetscFunctionReturn(PETSC_SUCCESS); 232 } 233 /*---------------------------------------------------------------------------------------------*/ 234 /*TEST 235 test: 236 suffix: 1 237 requires: !single 238 output_file: output/ex3span_1.out 239 args: -ts_monitor -ts_adapt_type none -restart 240 args: -dtpost 0.1127 -D 0.0015 -dir 0 -ts_max_time 9.8 -ts_dt 0.18 241 nsize: 1 242 243 test: 244 suffix: 1single 245 requires: single 246 output_file: output/ex3span_1single.out 247 args: -ts_monitor -ts_adapt_type none -restart -ts_event_dt_min 1e-6 248 args: -dtpost 0.1127 -D 0.0015 -dir 0 -ts_max_time 9.8 -ts_dt 0.18 249 nsize: 1 250 251 test: 252 suffix: 2 253 output_file: output/ex3span_2.out 254 args: -ts_event_dt_min 1e-6 -dtpost 1 -term 0 -ts_max_time 9.61 255 nsize: 1 256 257 test: 258 suffix: fin8 259 output_file: output/ex3span_fin8.out 260 args: -ts_max_time 8 261 args: -ts_monitor -ts_event_dt_min 1e-6 262 args: -ts_adapt_type {{none basic}} 263 args: -dtpost 0.1127 264 args: -D {{0.0015 0.03}} 265 args: -dir {{0 1 -1}} 266 args: -ts_dt 0.302 267 args: -ts_type {{beuler rk alpha rosw bdf}} 268 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 269 nsize: 1 270 271 test: 272 suffix: fin81 273 output_file: output/ex3span_fin8.1.out 274 args: -ts_max_time 8.1 275 args: -ts_monitor -ts_event_dt_min 1e-6 276 args: -ts_adapt_type {{none basic}} 277 args: -dtpost 0.1125 278 args: -D {{0.0015 0.03}} 279 args: -dir {{0 1 -1}} 280 args: -ts_dt 0.302 281 args: -ts_type {{beuler rk alpha rosw bdf}} 282 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 283 nsize: 2 284 285 test: 286 suffix: fin821 287 output_file: output/ex3span_fin8.21.out 288 args: -ts_max_time 8.21 289 args: -ts_monitor -ts_event_dt_min 1e-6 290 args: -ts_adapt_type {{none basic}} 291 args: -dtpost 0 292 args: -D {{0.0015 0.03}} 293 args: -dir {{0 1 -1}} 294 args: -ts_dt 0.302 295 args: -ts_type {{beuler rk alpha rosw bdf}} 296 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 297 nsize: 1 298 299 test: 300 suffix: fin899 301 output_file: output/ex3span_fin8.99.out 302 args: -ts_max_time 8.99 303 args: -ts_monitor -ts_event_dt_min 1e-6 304 args: -ts_adapt_type {{none basic}} 305 args: -dtpost 0.1127 306 args: -D {{0.0015 0.03}} 307 args: -dir {{0 1 -1}} 308 args: -ts_dt 0.302 309 args: -ts_type {{beuler rk alpha rosw bdf}} 310 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 311 nsize: 4 312 313 test: 314 suffix: fin9 315 output_file: output/ex3span_fin9.out 316 args: -ts_max_time 9 317 args: -ts_monitor -ts_event_dt_min 1e-6 318 args: -ts_adapt_type {{none basic}} 319 args: -dtpost 0 320 args: -D {{0.0015 0.03}} 321 args: -dir {{0 1 -1}} 322 args: -ts_dt 0.202 323 args: -ts_type {{beuler rk alpha rosw bdf}} 324 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 325 nsize: 2 326 327 test: 328 suffix: fin901 329 output_file: output/ex3span_fin9.01.out 330 args: -ts_max_time 9.01 331 args: -ts_monitor -ts_event_dt_min 1e-6 332 args: -ts_adapt_type {{none basic}} 333 args: -dtpost 0.1127 334 args: -D {{0.0015 0.03}} 335 args: -dir {{0 1 -1}} 336 args: -ts_dt 0.6 337 args: -ts_type {{beuler rk alpha rosw bdf}} 338 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 339 nsize: 1 340 341 test: 342 suffix: fin904 343 output_file: output/ex3span_fin9.04.out 344 args: -ts_max_time 9.04 345 args: -ts_monitor -ts_event_dt_min 1e-6 346 args: -ts_adapt_type {{none basic}} 347 args: -dtpost 0 348 args: -D {{0.0015 0.03}} 349 args: -dir {{0 1 -1}} 350 args: -ts_dt 0.102 351 args: -ts_type {{beuler rk alpha rosw bdf}} 352 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 353 nsize: 1 354 355 test: 356 suffix: fin905 357 output_file: output/ex3span_fin9.05.out 358 args: -ts_max_time 9.05 359 args: -ts_monitor -ts_event_dt_min 1e-6 360 args: -ts_adapt_type {{none basic}} 361 args: -dtpost 0.1127 362 args: -D {{0.0015 0.03}} 363 args: -dir {{0 1 -1}} 364 args: -ts_dt 0.402 365 args: -ts_type {{beuler rk alpha rosw bdf}} 366 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 367 nsize: 2 368 369 test: 370 suffix: fin905etc 371 output_file: output/ex3span_fin9.05.out 372 args: -ts_max_time {{9.06 9.07 9.1 9.21 9.5 9.99 10 11}} 373 args: -ts_monitor -ts_event_dt_min 1e-6 374 args: -ts_adapt_type {{none basic}} 375 args: -dtpost 0.11 376 args: -D 0.0025 377 args: -dir {{0 -1}} 378 args: -ts_dt 0.3025 379 args: -ts_type {{beuler rk alpha rosw bdf}} 380 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 381 nsize: 2 382 383 test: 384 suffix: fin906 385 output_file: output/ex3span_fin9.06.out 386 args: -ts_max_time 9.06 387 args: -ts_monitor -ts_event_dt_min 1e-6 388 args: -ts_adapt_type {{none basic}} 389 args: -dtpost 0.1121 390 args: -D {{0.0015 0.02}} 391 args: -dir 1 392 args: -ts_dt 0.302 393 args: -ts_type {{beuler rk alpha rosw bdf}} 394 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 395 nsize: 1 396 397 test: 398 suffix: fin91 399 output_file: output/ex3span_fin9.1.out 400 args: -ts_max_time 9.1 401 args: -ts_monitor -ts_event_dt_min 1e-6 402 args: -ts_adapt_type {{none basic}} 403 args: -dtpost 0.09 404 args: -D {{0.0015 0.03}} 405 args: -dir 1 406 args: -ts_dt 0.302 407 args: -ts_type {{beuler rk alpha rosw bdf}} 408 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 409 nsize: 2 410 411 test: 412 suffix: fin921 413 output_file: output/ex3span_fin9.21.out 414 args: -ts_max_time 9.21 415 args: -ts_monitor -ts_event_dt_min 1e-6 416 args: -ts_adapt_type {{none basic}} 417 args: -dtpost 0.7 418 args: -D {{0.0015 0.02}} 419 args: -dir 1 420 args: -ts_dt 0.22 421 args: -ts_type {{beuler rk alpha rosw bdf}} 422 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 423 nsize: 1 424 425 test: 426 suffix: fin95 427 output_file: output/ex3span_fin9.5.out 428 args: -ts_max_time 9.5 429 args: -ts_monitor -ts_event_dt_min 1e-6 430 args: -ts_adapt_type {{none basic}} 431 args: -dtpost 0 432 args: -D {{0.0015 0.0135}} 433 args: -dir 1 434 args: -ts_dt 0.502 435 args: -ts_type {{beuler rk alpha rosw bdf}} 436 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 437 nsize: 4 438 439 test: 440 suffix: fin10 441 output_file: output/ex3span_fin10.out 442 args: -ts_max_time 10 443 args: -ts_monitor -ts_event_dt_min 1e-6 444 args: -ts_adapt_type {{none basic}} 445 args: -dtpost 0.17 446 args: -D {{0.0015 0.02}} 447 args: -dir 1 448 args: -ts_dt 0.422 449 args: -ts_type {{beuler rk alpha rosw bdf}} 450 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 451 nsize: 2 452 453 test: 454 suffix: fin11 455 output_file: output/ex3span_fin11.out 456 args: -ts_max_time 11 457 args: -ts_monitor -ts_event_dt_min 1e-6 458 args: -ts_adapt_type {{none basic}} 459 args: -dtpost 0 460 args: -D {{0.0015 0.0016}} 461 args: -dir 1 462 args: -ts_dt 0.2718281828 463 args: -ts_type {{beuler rk alpha rosw bdf}} 464 filter: grep "TS dt" | tail -n 1 | cut -f6 -d " " 465 nsize: 4 466 467 TEST*/ 468