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 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 (default), 1, -1\n" 17 "-flg : additional output in Postevent (default: nothing)\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 (default: no)\n" 20 "-term : flag to terminate at 9.05 event (true by default)\n" 21 "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n" 22 " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n" 23 " if x == 0, nothing happens (default)\n" 24 "-D z : a small real number to define additional TimeSpan points (default = 0.02)\n" 25 "-dt2_at6 t : second time step set after event at t=6 (if nothing is specified, no action is done)\n" 26 "-mult7 m : after event at t=7, the linear system coeffs '0.2' are multiplied by m (default = 1.0)\n"; 27 28 #define MAX_NFUNC 100 // max event functions per rank 29 #define MAX_NEV 5000 // max zero crossings for each rank 30 31 typedef struct { 32 PetscMPIInt rank, size; 33 PetscReal pi; 34 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 35 PetscReal evres[MAX_NEV]; // times of found zero-crossings 36 PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking 37 PetscInt cnt; // counter 38 PetscInt cntref; // actual length of 'ref' on the given rank 39 PetscBool flg; // flag for additional print in PostEvent 40 PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default) 41 PetscBool restart; // flag for TSRestartStep() in PostEvent 42 PetscBool term; // flag to terminate at 9.05 event 43 PetscReal dtpost; // first post-event step 44 PetscReal dt2_at6; // second time step set after event at t=6 45 PetscReal mult7; // multiplier for coeffs at t=7 46 PetscInt postcnt; // counter for PostEvent calls 47 Mat A; // system matrix 48 PetscInt m; // local size of A 49 } AppCtx; 50 51 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 52 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 53 PetscErrorCode Fill_mat(PetscReal coeff, PetscInt m, Mat A); // Fills the system matrix (2*2) 54 55 int main(int argc, char **argv) 56 { 57 TS ts; 58 Vec sol; 59 PetscInt n, dir0; 60 PetscReal tol = 1e-7, D = 0.02; 61 PetscInt dir[MAX_NFUNC]; 62 PetscBool term[MAX_NFUNC], match; 63 PetscScalar *x; 64 PetscReal tspan[28], dtlast, tlast, tlast_expected, maxtime; 65 AppCtx ctx; 66 TSConvergedReason reason; 67 TSAdapt adapt; 68 69 PetscFunctionBeginUser; 70 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 71 setbuf(stdout, NULL); 72 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 73 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 74 ctx.pi = PetscAcosReal(-1.0); 75 ctx.cnt = 0; 76 ctx.cntref = 0; 77 ctx.flg = PETSC_FALSE; 78 ctx.errtol = 1e-5; 79 ctx.restart = PETSC_FALSE; 80 ctx.term = PETSC_TRUE; 81 ctx.dtpost = 0; 82 ctx.dt2_at6 = -2; 83 ctx.mult7 = 1.0; 84 ctx.postcnt = 0; 85 ctx.m = 0; 86 87 // The linear problem has a 2*2 matrix. The matrix is constant 88 if (ctx.rank == 0) ctx.m = 2; 89 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, ctx.m, ctx.m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &ctx.A)); 90 PetscCallBack("Fill_mat", Fill_mat(0.2, ctx.m, ctx.A)); 91 PetscCall(MatCreateVecs(ctx.A, &sol, NULL)); 92 PetscCall(VecGetArray(sol, &x)); 93 if (ctx.rank == 0) { // initial conditions 94 x[0] = 0; // sin(0) 95 x[1] = 1; // cos(0) 96 } 97 PetscCall(VecRestoreArray(sol, &x)); 98 99 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 100 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 101 102 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 103 PetscCall(TSSetRHSJacobian(ts, ctx.A, ctx.A, TSComputeRHSJacobianConstant, NULL)); 104 105 PetscCall(TSSetTimeStep(ts, 0.099)); 106 PetscCall(TSSetType(ts, TSBEULER)); 107 PetscCall(TSSetMaxSteps(ts, 10000)); 108 PetscCall(TSSetMaxTime(ts, 10.0)); 109 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 110 111 // Set the event handling 112 dir0 = 0; 113 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction 114 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 115 PetscCall(PetscOptionsGetReal(NULL, NULL, "-errtol", &ctx.errtol, NULL)); // error tolerance for located events 116 PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep() 117 PetscCall(PetscOptionsGetBool(NULL, NULL, "-term", &ctx.term, NULL)); // flag to terminate at 9.05 event 118 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step 119 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dt2_at6", &ctx.dt2_at6, NULL)); // second time step set after event at t=6 120 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mult7", &ctx.mult7, NULL)); // multiplier for coeffs at t=7 121 PetscCall(PetscOptionsGetReal(NULL, NULL, "-D", &D, NULL)); // small number for tspan 122 123 n = 0; // event counter 124 if (ctx.rank == 0) { // first event -- on rank-0 125 dir[n] = dir0; 126 term[n++] = PETSC_FALSE; 127 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.05; 128 } 129 if (ctx.rank == ctx.size - 1) { // second event (with optional termination) -- on last rank 130 dir[n] = dir0; 131 term[n++] = ctx.term; 132 if (dir0 <= 0) ctx.ref[ctx.cntref++] = 9.05; 133 } 134 if (ctx.rank == 1 % ctx.size) { // third event -- on rank = 1%ctx.size 135 dir[n] = dir0; 136 term[n++] = PETSC_FALSE; 137 138 for (PetscInt i = 1; i < MAX_NEV - 2; i++) { 139 if (i % 2 == 1 && dir0 <= 0) ctx.ref[ctx.cntref++] = i; 140 if (i % 2 == 0 && dir0 >= 0) ctx.ref[ctx.cntref++] = i; 141 } 142 } 143 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref)); 144 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 145 PetscCall(TSSetEventTolerances(ts, tol, NULL)); 146 147 // Set the time span 148 for (PetscInt i = 0; i < 10; i++) { 149 tspan[2 * i] = 0.01 + i + (i == 7 ? -0.02 : 0); 150 tspan[2 * i + 1] = 0.21 + i; 151 } 152 tspan[20] = 3; 153 tspan[21] = 4; 154 tspan[22] = 4 + D; 155 tspan[23] = 5 - D; 156 tspan[24] = 5; 157 tspan[25] = 6 - D; 158 tspan[26] = 6; 159 tspan[27] = 6 + D; 160 PetscCall(PetscSortReal(28, tspan)); 161 PetscCall(TSSetTimeSpan(ts, 28, tspan)); 162 PetscCall(TSSetFromOptions(ts)); 163 164 // Solution 165 PetscCall(TSSolve(ts, sol)); 166 PetscCall(TSGetConvergedReason(ts, &reason)); 167 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CONVERGED REASON: %" PetscInt_FMT " (TS_CONVERGED_EVENT == %" PetscInt_FMT ")\n", (PetscInt)reason, (PetscInt)TS_CONVERGED_EVENT)); 168 169 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"] 170 for (PetscInt j = 0; j < ctx.cnt; j++) { 171 PetscReal err = 10.0; 172 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]); 173 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")); 174 } 175 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 176 177 // print the final time and step 178 PetscCall(TSGetTime(ts, &tlast)); 179 PetscCall(TSGetTimeStep(ts, &dtlast)); 180 PetscCall(TSGetAdapt(ts, &adapt)); 181 PetscCall(PetscObjectTypeCompare((PetscObject)adapt, TSADAPTNONE, &match)); 182 183 PetscCall(TSGetMaxTime(ts, &maxtime)); 184 tlast_expected = ((dir0 == 1 || !ctx.term) ? maxtime : PetscMin(maxtime, 9.05)); 185 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final time = %g, max time = %g, %s\n", (double)tlast, (double)maxtime, PetscAbsReal(tlast - tlast_expected) < ctx.errtol ? "pass" : "fail")); 186 187 if (match) { 188 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Adapt = none\n")); 189 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Last dt = %g\n", (double)dtlast)); 190 } 191 192 PetscCall(MatDestroy(&ctx.A)); 193 PetscCall(TSDestroy(&ts)); 194 PetscCall(VecDestroy(&sol)); 195 196 PetscCall(PetscFinalize()); 197 return 0; 198 } 199 200 /* 201 User callback for defining the event-functions 202 */ 203 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 204 { 205 PetscInt n = 0; 206 AppCtx *Ctx = (AppCtx *)ctx; 207 208 PetscFunctionBeginUser; 209 // for the test purposes, event-functions are defined based on t 210 // first event -- on rank-0 211 if (Ctx->rank == 0) { 212 if (t < 2.05) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.05, 12)); 213 else gval[n++] = 0.5; 214 } 215 216 // second event -- on last rank 217 if (Ctx->rank == Ctx->size - 1) { 218 if (t > 8.05) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.05, 12)); 219 else gval[n++] = 0.25; 220 } 221 222 // third event -- on rank = 1%ctx.size 223 if (Ctx->rank == 1 % Ctx->size) { gval[n++] = PetscSinReal(Ctx->pi * t); } 224 PetscFunctionReturn(PETSC_SUCCESS); 225 } 226 227 /* 228 User callback for the post-event stuff 229 */ 230 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 231 { 232 AppCtx *Ctx = (AppCtx *)ctx; 233 PetscBool mat_changed = PETSC_FALSE; 234 235 PetscFunctionBeginUser; 236 if (Ctx->flg) { 237 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 238 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 239 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 240 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 241 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 242 } 243 244 if (Ctx->cnt + nev_zero < MAX_NEV) 245 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing 246 247 #ifdef NEW_VERSION 248 Ctx->postcnt++; // sync 249 if (Ctx->dtpost > 0) { 250 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 251 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE)); 252 } 253 #endif 254 255 // t==6: set the second post-event step 256 if (PetscAbsReal(t - 6.0) < 0.01 && Ctx->dt2_at6 != -2) PetscCall(TSSetPostEventSecondStep(ts, Ctx->dt2_at6)); 257 258 // t==7: change the system matrix 259 if (PetscAbsReal(t - 7.0) < 0.01 && Ctx->mult7 != 1) { 260 PetscCallBack("Fill_mat", Fill_mat(0.2 * Ctx->mult7, Ctx->m, Ctx->A)); 261 PetscCall(TSSetRHSJacobian(ts, Ctx->A, Ctx->A, TSComputeRHSJacobianConstant, NULL)); 262 mat_changed = PETSC_TRUE; 263 } 264 265 if (Ctx->restart || mat_changed) PetscCall(TSRestartStep(ts)); 266 PetscFunctionReturn(PETSC_SUCCESS); 267 } 268 269 /* 270 Fills the system matrix (2*2) 271 */ 272 PetscErrorCode Fill_mat(PetscReal coeff, PetscInt m, Mat A) 273 { 274 PetscInt inds[2]; 275 PetscScalar vals[4]; 276 277 PetscFunctionBeginUser; 278 inds[0] = 0; 279 inds[1] = 1; 280 vals[0] = 0; 281 vals[1] = coeff; 282 vals[2] = -coeff; 283 vals[3] = 0; 284 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 285 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 286 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 287 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 288 PetscFunctionReturn(PETSC_SUCCESS); 289 } 290 /*---------------------------------------------------------------------------------------------*/ 291 /* 292 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1, 293 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to 294 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed, 295 simply remove this option altogether. This will result in using the defaults 296 (which is PETSC_DECIDE). 297 */ 298 /*TEST 299 test: 300 suffix: 1 301 requires: !single 302 output_file: output/ex3span_1.out 303 args: -ts_monitor -ts_adapt_type none -restart 304 args: -dtpost 0.1127 -D 0.0015 -dir 0 -ts_max_time 9.8 -ts_dt 0.18 305 nsize: 1 306 307 test: 308 suffix: 1single 309 requires: single 310 output_file: output/ex3span_1single.out 311 args: -ts_monitor -ts_adapt_type none -restart -ts_event_dt_min 1e-6 312 args: -dtpost 0.1127 -D 0.0015 -dir 0 -ts_max_time 9.8 -ts_dt 0.18 313 nsize: 1 314 315 test: 316 suffix: 2 317 output_file: output/ex3span_2.out 318 args: -ts_event_dt_min 1e-6 -dtpost 1 -term 0 -ts_max_time 9.61 319 nsize: 1 320 321 test: 322 suffix: 3none 323 output_file: output/ex3span_3none.out 324 args: -ts_event_dt_min 1e-6 -ts_adapt_type none -dir 0 325 args: -ts_event_post_event_step {{-1 0.11}} 326 args: -ts_event_post_event_second_step 0.12 327 args: -dt2_at6 {{-2 0.08 0.15}} 328 nsize: 3 329 330 test: 331 suffix: 3basic 332 output_file: output/ex3span_3basic.out 333 args: -ts_event_dt_min 1e-6 -ts_adapt_type basic -dir 0 334 args: -ts_event_post_event_step {{-1 0.11}} 335 args: -ts_event_post_event_second_step 0.12 336 args: -dt2_at6 {{-2 0.08 0.15}} 337 args: -mult7 {{1 2}} 338 nsize: 2 339 340 test: 341 suffix: fin 342 output_file: output/ex3span_fin.out 343 args: -ts_max_time {{8.21 8.99 9 9.04 9.05 9.06 9.21 9.99 12}} 344 args: -ts_event_dt_min 1e-6 345 args: -ts_adapt_type {{none basic}} 346 args: -dtpost 0.1125 347 args: -D 0.0025 348 args: -dir {{0 -1 1}} 349 args: -ts_dt 0.3025 350 args: -ts_type {{rk bdf}} 351 filter: grep "Final time =" 352 filter_output: grep "Final time =" 353 nsize: 2 354 355 test: 356 suffix: adaptmonitor 357 requires: !single 358 output_file: output/ex3span_adaptmonitor.out 359 args: -ts_adapt_monitor -dir 1 360 nsize: 1 361 TEST*/ 362