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 "Using one or several event functions (on rank-0)\n" 10 "This program is mostly intended to test the Anderson-Bjorck iteration with challenging event-functions\n" 11 "Options:\n" 12 "-dir d : zero-crossing direction for events\n" 13 "-flg : additional output in Postevent\n" 14 "-errtol e : error tolerance, for printing 'pass/fail' for located events (1e-5 by default)\n" 15 "-restart : flag for TSRestartStep() in PostEvent\n" 16 "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n" 17 " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n" 18 " if x == 0, nothing happens\n" 19 "-func F : selects the event function [0, ..., 11], if F == -1 (default) is set, all event functions are taken\n"; 20 21 #define MAX_NFUNC 100 // max event functions per rank 22 #define MAX_NEV 5000 // max zero crossings for each rank 23 24 typedef struct { 25 PetscMPIInt rank, size; 26 PetscReal pi; 27 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 28 PetscReal evres[MAX_NEV]; // times of found zero-crossings 29 PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking 30 PetscInt cnt; // counter 31 PetscInt cntref; // actual length of 'ref' on the given rank 32 PetscBool flg; // flag for additional print in PostEvent 33 PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default) 34 PetscBool restart; // flag for TSRestartStep() in PostEvent 35 PetscReal dtpost; // post-event step 36 PetscInt postcnt; // counter for PostEvent calls 37 PetscInt F; // event-function index 38 PetscInt Fnum; // total available event functions 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 PetscInt dir[MAX_NFUNC], inds[2]; 51 PetscBool term[MAX_NFUNC]; 52 PetscBool pass = PETSC_TRUE; 53 PetscScalar *x, vals[4]; 54 AppCtx ctx; 55 56 PetscFunctionBeginUser; 57 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 58 setbuf(stdout, NULL); 59 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 60 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 61 ctx.pi = PetscAcosReal(-1.0); 62 ctx.cnt = 0; 63 ctx.cntref = 0; 64 ctx.flg = PETSC_FALSE; 65 ctx.errtol = 1e-5; 66 ctx.restart = PETSC_FALSE; 67 ctx.dtpost = 0; 68 ctx.postcnt = 0; 69 ctx.F = -1; 70 ctx.Fnum = 12; 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(TSSetTime(ts, 0.03)); 101 PetscCall(TSSetTimeStep(ts, 0.1)); 102 PetscCall(TSSetType(ts, TSBEULER)); 103 PetscCall(TSSetMaxSteps(ts, 10000)); 104 PetscCall(TSSetMaxTime(ts, 4.0)); 105 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 106 PetscCall(TSSetFromOptions(ts)); 107 108 // Set the event handling 109 dir0 = 0; 110 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction 111 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 112 PetscCall(PetscOptionsGetReal(NULL, NULL, "-errtol", &ctx.errtol, NULL)); // error tolerance for located events 113 PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep() 114 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step 115 PetscCall(PetscOptionsGetInt(NULL, NULL, "-F", &ctx.F, NULL)); // event-function index 116 PetscCheck(ctx.F >= -1 && ctx.F < ctx.Fnum, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Value of 'F' is out of range"); 117 118 n = 0; // event counter 119 if (ctx.rank == 0) { // all events -- on rank-0 120 if (ctx.F == -1) 121 for (n = 0; n < ctx.Fnum; n++) { // all event-functions 122 dir[n] = dir0; 123 term[n] = PETSC_FALSE; 124 } 125 else { // single event-function 126 dir[n] = dir0; 127 term[n++] = PETSC_FALSE; 128 } 129 130 // set the reference values 131 if (ctx.F == 0 || ctx.F == -1) { 132 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 2.0; 133 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 4.0; 134 if (dir0 <= 0) ctx.ref[ctx.cntref++] = 1.0; 135 if (dir0 <= 0) ctx.ref[ctx.cntref++] = 3.0; 136 } 137 if (ctx.F == 1 || ctx.F == -1) 138 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.0; 139 if (ctx.F == 2 || ctx.F == -1) 140 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.0; 141 if (ctx.F == 3 || ctx.F == -1) 142 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.696812386809752; 143 if (ctx.F == 4 || ctx.F == -1) 144 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.1; 145 if (ctx.F == 5 || ctx.F == -1) 146 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 1.02; 147 if (ctx.F == 6 || ctx.F == -1) 148 if (dir0 <= 0) ctx.ref[ctx.cntref++] = 0.90478821787302; 149 if (ctx.F == 7 || ctx.F == -1) 150 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.05; 151 if (ctx.F == 8 || ctx.F == -1) 152 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.552704666678489; 153 if (ctx.F == 9 || ctx.F == -1) 154 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.399422291710969; 155 if (ctx.F == 10 || ctx.F == -1) { 156 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 0.874511220376091; 157 if (dir0 <= 0) ctx.ref[ctx.cntref++] = 2.388359335869107; 158 } 159 if (ctx.F == 11 || ctx.F == -1) 160 if (dir0 >= 0) ctx.ref[ctx.cntref++] = 2.5; 161 } 162 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref)); 163 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 164 165 // Solution 166 PetscCall(TSSolve(ts, sol)); 167 168 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"] 169 for (PetscInt j = 0; j < ctx.cnt; j++) { 170 PetscReal err = 10.0; 171 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]); 172 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")); 173 pass = (pass && err < ctx.errtol ? PETSC_TRUE : PETSC_FALSE); 174 } 175 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 176 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "This test: %s\n", pass ? "PASSED" : "FAILED")); 177 178 PetscCall(MatDestroy(&A)); 179 PetscCall(TSDestroy(&ts)); 180 PetscCall(VecDestroy(&sol)); 181 182 PetscCall(PetscFinalize()); 183 return 0; 184 } 185 186 /* 187 User callback for defining the event-functions 188 */ 189 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 190 { 191 PetscInt n = 0; 192 AppCtx *Ctx = (AppCtx *)ctx; 193 194 PetscFunctionBeginUser; 195 // for the test purposes, event-functions are defined based on t 196 // all events -- on rank-0 197 if (Ctx->rank == 0) { 198 if (Ctx->F == 0 || Ctx->F == -1) gval[n++] = PetscSinReal(Ctx->pi * t) / Ctx->pi; // FUNC-0, roots 1, 2, 3, 4 199 if (Ctx->F == 1 || Ctx->F == -1) gval[n++] = PetscLogReal(t); // FUNC-2, root 1 200 if (Ctx->F == 2 || Ctx->F == -1) { // FUNC-3, root 1 201 if (t < 2) gval[n++] = (1 - PetscPowReal(t - 2, 12)) / 12.0; 202 else gval[n++] = 1 / 12.0; 203 } 204 if (Ctx->F == 3 || Ctx->F == -1) gval[n++] = t - PetscExpReal(PetscSinReal(t)) + 1; // FUNC-5, root 1.69681 205 if (Ctx->F == 4 || Ctx->F == -1) gval[n++] = (1e10 * PetscPowReal(t, 1 / t) - 1) / 100; // FUNC-6, root 0.1 206 if (Ctx->F == 5 || Ctx->F == -1) gval[n++] = PetscLogReal(t - 0.02) * PetscLogReal(t - 0.02) * PetscSignReal(t - 1.02) * 1e7; // FUNC-7, root 1.02 207 if (Ctx->F == 6 || Ctx->F == -1) gval[n++] = 4 * PetscCosReal(t) - PetscExpReal(t); // FUNC-14, root 0.904788 208 if (Ctx->F == 7 || Ctx->F == -1) gval[n++] = (20.0 * t - 1) / (19.0 * t) / 10; // FUNC-15, root 0.05 209 if (Ctx->F == 8 || Ctx->F == -1) gval[n++] = ((t - 1) * PetscExpReal(-20 * t) + PetscPowReal(t, 20)) * 1e4; // FUNC-16, root 0.5527 210 if (Ctx->F == 9 || Ctx->F == -1) gval[n++] = (t * t * (t * t / 3.0 + PetscSqrtReal(2.0) * PetscSinReal(t)) - PetscSqrtReal(3.0) / 18) * 10; // FUNC-17, root 0.399 211 if (Ctx->F == 10 || Ctx->F == -1) gval[n++] = ((t * t + 1) * PetscSinReal(t) - PetscExpReal(PetscSqrtReal(t)) * (t - 1) * (t * t - 5)) / 10; // FUNC-18, roots 0.87, 2.388 212 if (Ctx->F == 11 || Ctx->F == -1) gval[n++] = 2 * t - 5; // FUNC-21, root 2.5 213 } 214 PetscFunctionReturn(PETSC_SUCCESS); 215 } 216 217 /* 218 User callback for the post-event stuff 219 */ 220 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 221 { 222 AppCtx *Ctx = (AppCtx *)ctx; 223 224 PetscFunctionBeginUser; 225 if (Ctx->flg) { 226 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 227 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 228 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 229 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 230 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 231 } 232 233 if (Ctx->cnt + nev_zero < MAX_NEV) 234 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing 235 236 #ifdef NEW_VERSION 237 Ctx->postcnt++; // sync 238 if (Ctx->dtpost > 0) { 239 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 240 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE)); 241 } 242 #endif 243 244 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 245 PetscFunctionReturn(PETSC_SUCCESS); 246 } 247 /*---------------------------------------------------------------------------------------------*/ 248 /* 249 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1, 250 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to 251 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed, 252 simply remove this option altogether. This will result in using the defaults 253 (which is PETSC_DECIDE). 254 */ 255 /*TEST 256 test: 257 suffix: 0 258 requires: !single 259 output_file: output/ex4_0.out 260 args: -dir 0 261 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-10 262 args: -ts_dt 0.25 263 args: -restart 0 264 args: -ts_event_tol {{1e-8 1e-15}} 265 args: -errtol 1e-7 266 args: -ts_adapt_type {{none basic}} 267 args: -dtpost 0 268 args: -ts_event_post_event_step -1 269 args: -ts_type rk 270 nsize: 2 271 filter: sort 272 filter_output: sort 273 274 test: 275 suffix: 0single 276 requires: single 277 output_file: output/ex4_0single.out 278 args: -dir 0 279 args: -ts_adapt_dt_min 1e-6 -ts_event_dt_min 1e-6 280 args: -ts_dt 0.3 281 args: -ts_event_tol {{1e-7 1e-10}} 282 args: -ts_adapt_type {{none basic}} 283 args: -dtpost 0.23 284 args: -ts_event_post_event_step -1 285 args: -ts_type beuler 286 nsize: 3 287 filter: sort 288 filter_output: sort 289 290 test: 291 suffix: F7 292 output_file: output/ex4_F7.out 293 args: -dir 0 294 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6 295 args: -ts_dt 0.4 296 args: -F 7 297 args: -ts_event_tol {{1e-8 1e-15}} 298 args: -ts_adapt_type {{none basic}} 299 args: -ts_type rk 300 nsize: 1 301 302 test: 303 suffix: F7revisit 304 output_file: output/ex4_F7revisit.out 305 args: -ts_event_monitor -F 7 -ts_dt 0.04 -ts_event_dt_min 0.016 -errtol 0.005 306 nsize: 1 307 308 test: 309 suffix: 2all 310 output_file: output/ex4_2.out 311 args: -dir 0 312 args: -F {{-1 0 1 2 3 4 5 6 7 8 9 10 11}} 313 args: -ts_event_dt_min 1e-6 -ts_dt 0.4 -ts_event_tol 1e-8 314 args: -ts_adapt_type {{none basic}} 315 args: -dtpost 0.35 316 args: -ts_type rk 317 filter: grep "This test" 318 nsize: 1 319 320 test: 321 suffix: 2pos 322 output_file: output/ex4_2.out 323 args: -dir 1 324 args: -F {{-1 0 1 2 3 4 5 7 8 9 10 11}} 325 args: -ts_event_dt_min 1e-6 -ts_dt 0.4 -ts_event_tol 1e-8 326 args: -ts_adapt_type none 327 args: -dtpost 0.34 328 args: -ts_type beuler 329 filter: grep "This test" 330 nsize: 1 331 332 test: 333 suffix: 2neg 334 output_file: output/ex4_2.out 335 args: -dir -1 336 args: -F {{-1 0 6 10}} 337 args: -ts_event_dt_min 1e-6 -ts_dt 0.4 -ts_event_tol 1e-8 338 args: -ts_adapt_type {{none basic}} 339 args: -dtpost 0.33 340 args: -ts_type rk 341 filter: grep "This test" 342 nsize: 1 343 344 test: 345 suffix: pos 346 output_file: output/ex4_pos.out 347 args: -dir 1 348 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6 349 args: -ts_dt 0.4 350 args: -restart 0 351 args: -ts_event_tol {{1e-8 1e-15}} 352 args: -ts_adapt_type {{none basic}} 353 args: -dtpost 0.25 354 args: -ts_event_post_event_step -1 355 args: -ts_type {{beuler rk}} 356 nsize: 1 357 358 test: 359 suffix: neg 360 output_file: output/ex4_neg.out 361 args: -dir -1 362 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6 363 args: -ts_dt 0.4 364 args: -restart 1 365 args: -ts_event_tol {{1e-8 1e-15}} 366 args: -ts_adapt_type {{none basic}} 367 args: -dtpost 0.25 368 args: -ts_event_post_event_step 0.35 369 args: -ts_type rk 370 nsize: 2 371 filter: sort 372 filter_output: sort 373 TEST*/ 374