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 16 event functions:\n" 10 "7 polynomials (dir=+1) with zeros: 1+2^i, i=-3,...3, on ranks=(i+3)%size\n" 11 "7 polynomials (dir=-1) with zeros: 1+(8-2^i), i=-3,...3, on ranks=(i+3)%size\n" 12 "(t-5)^2 * sin(pi*t), with zeros = 1,2,...10, on rank-0\n" 13 " 0.5 * cos(pi*t), with zeros = 0.5,1.5,...9.5, on last-rank\n" 14 "Options:\n" 15 "-dir d : zero-crossing direction for events\n" 16 "-flg : additional output in Postevent\n" 17 "-errtol e : error tolerance, for printing 'pass/fail' for located events (1e-5 by default)\n" 18 "-restart : flag for TSRestartStep() in PostEvent\n" 19 "-dtpost x : if x > 0, then on even PostEvent calls 1st-post-event-step = x is set,\n" 20 " on odd PostEvent calls 1st-post-event-step = PETSC_DECIDE is set,\n" 21 " if x == 0, nothing happens\n"; 22 23 #define MAX_NFUNC 100 // max event functions per rank 24 #define MAX_NEV 5000 // max zero crossings for each rank 25 26 typedef struct { 27 PetscMPIInt rank, size; 28 PetscReal pi; 29 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 30 PetscReal evres[MAX_NEV]; // times of found zero-crossings 31 PetscReal ref[MAX_NEV]; // reference times of zero-crossings, for checking 32 PetscInt cnt; // counter 33 PetscInt cntref; // actual length of 'ref' on the given rank 34 PetscBool flg; // flag for additional print in PostEvent 35 PetscReal errtol; // error tolerance, for printing 'pass/fail' for located events (1e-5 by default) 36 PetscBool restart; // flag for TSRestartStep() in PostEvent 37 PetscReal dtpost; // post-event step 38 PetscInt postcnt; // counter for PostEvent calls 39 PetscReal vtol[MAX_NFUNC]; // vtol array, with extra storage 40 PetscInt dir0; // desired zero-crossing direction 41 } AppCtx; 42 43 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 44 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 45 static inline void SetVtols(PetscMPIInt rank, PetscMPIInt size, PetscReal tol0, PetscReal tolsin, PetscReal *vtol); // helper function to fill vtol[] 46 47 int main(int argc, char **argv) 48 { 49 TS ts; 50 Mat A; 51 Vec sol; 52 PetscInt n, m = 0; 53 PetscInt dir[MAX_NFUNC], inds[2]; 54 PetscBool term[MAX_NFUNC]; 55 PetscScalar *x, vals[4]; 56 PetscReal aux; 57 AppCtx ctx; 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.cntref = 0; 67 ctx.flg = PETSC_FALSE; 68 ctx.errtol = 1e-5; 69 ctx.restart = PETSC_FALSE; 70 ctx.dtpost = 0; 71 ctx.postcnt = 0; 72 73 // The linear problem has a 2*2 matrix. The matrix is constant 74 if (ctx.rank == 0) m = 2; 75 inds[0] = 0; 76 inds[1] = 1; 77 vals[0] = 0; 78 vals[1] = 0.2; 79 vals[2] = -0.2; 80 vals[3] = 0; 81 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A)); 82 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 83 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 84 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 85 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 86 87 PetscCall(MatCreateVecs(A, &sol, NULL)); 88 PetscCall(VecGetArray(sol, &x)); 89 if (ctx.rank == 0) { // initial conditions 90 x[0] = 0; // sin(0) 91 x[1] = 1; // cos(0) 92 } 93 PetscCall(VecRestoreArray(sol, &x)); 94 95 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 96 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 97 98 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 99 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL)); 100 101 PetscCall(TSSetTimeStep(ts, 0.1)); 102 PetscCall(TSSetType(ts, TSBEULER)); 103 PetscCall(TSSetMaxSteps(ts, 10000)); 104 PetscCall(TSSetMaxTime(ts, 10.0)); 105 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 106 PetscCall(TSSetFromOptions(ts)); 107 108 // Set the event handling 109 ctx.dir0 = 0; 110 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &ctx.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 116 n = 0; // event counter 117 aux = 1.0 / 8.0; 118 for (PetscInt i = -3; i <= 3; i++) { // pos-polynomials 119 if (ctx.rank == (i + 3) % ctx.size) { 120 dir[n] = ctx.dir0; 121 term[n++] = PETSC_FALSE; 122 if (ctx.dir0 >= 0) ctx.ref[ctx.cntref++] = 1.0 + aux; 123 } 124 aux *= 2; 125 } 126 aux = 1.0 / 8.0; 127 for (PetscInt i = -3; i <= 3; i++) { // neg-polynomials 128 if (ctx.rank == (i + 3) % ctx.size) { 129 dir[n] = ctx.dir0; 130 term[n++] = PETSC_FALSE; 131 if (ctx.dir0 <= 0) ctx.ref[ctx.cntref++] = 9.0 - aux; 132 } 133 aux *= 2; 134 } 135 if (ctx.rank == 0) { // sin-event -- on rank-0 136 dir[n] = ctx.dir0; 137 term[n++] = PETSC_FALSE; 138 for (PetscInt i = 1; i < MAX_NEV / 2 - 10; i++) { 139 if (i % 2 == 1 && ctx.dir0 <= 0) ctx.ref[ctx.cntref++] = i; 140 if (i % 2 == 0 && ctx.dir0 >= 0) ctx.ref[ctx.cntref++] = i; 141 } 142 } 143 if (ctx.rank == ctx.size - 1) { // cos-event -- on last rank 144 dir[n] = ctx.dir0; 145 term[n++] = PETSC_FALSE; 146 for (PetscInt i = 1; i < MAX_NEV / 2 - 10; i++) { 147 if (i % 2 == 1 && ctx.dir0 <= 0) ctx.ref[ctx.cntref++] = i - 0.5; 148 if (i % 2 == 0 && ctx.dir0 >= 0) ctx.ref[ctx.cntref++] = i - 0.5; 149 } 150 } 151 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref)); 152 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 153 SetVtols(ctx.rank, ctx.size, 1e-8, 1e-8, ctx.vtol); 154 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, ctx.vtol)); 155 156 // Solution 157 PetscCall(TSSolve(ts, sol)); 158 159 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"] 160 for (PetscInt j = 0; j < ctx.cnt; j++) { 161 PetscReal err = 10.0; 162 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]); 163 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")); 164 } 165 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 166 167 PetscCall(MatDestroy(&A)); 168 PetscCall(TSDestroy(&ts)); 169 PetscCall(VecDestroy(&sol)); 170 171 PetscCall(PetscFinalize()); 172 return 0; 173 } 174 175 /* 176 User callback for defining the event-functions 177 */ 178 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 179 { 180 PetscInt n = 0; 181 AppCtx *Ctx = (AppCtx *)ctx; 182 PetscReal P; 183 184 PetscFunctionBeginUser; 185 // for the test purposes, event-functions are defined based on t 186 for (PetscInt i = -3; i <= 3; i++) { // pos-polynomials 187 if (Ctx->rank == (i + 3) % Ctx->size) { 188 P = PetscPowReal(2.0, i); 189 if (t < 2 + P) gval[n++] = 1 - PetscPowReal(2 + P - t, i + 5); 190 else gval[n++] = 1; 191 } 192 } 193 for (PetscInt i = -3; i <= 3; i++) { // neg-polynomials 194 if (Ctx->rank == (i + 3) % Ctx->size) { 195 P = PetscPowReal(2.0, i); 196 if (t > 8 - P) gval[n++] = 1 - PetscPowReal(t - 8 + P, i + 5); 197 else gval[n++] = 1; 198 } 199 } 200 if (Ctx->rank == 0) gval[n++] = (t - 5) * (t - 5) * PetscSinReal(Ctx->pi * t); // sin-event -- on rank-0 201 if (Ctx->rank == Ctx->size - 1) gval[n++] = 0.5 * PetscCosReal(Ctx->pi * t); // cos-event -- on last rank 202 PetscFunctionReturn(PETSC_SUCCESS); 203 } 204 205 /* 206 User callback for the post-event stuff 207 */ 208 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 209 { 210 AppCtx *Ctx = (AppCtx *)ctx; 211 212 PetscFunctionBeginUser; 213 if (Ctx->flg) { 214 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 215 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 216 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 217 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 218 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 219 } 220 221 if (Ctx->cnt + nev_zero < MAX_NEV) 222 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing 223 224 #ifdef NEW_VERSION 225 Ctx->postcnt++; // sync 226 if (Ctx->dtpost > 0) { 227 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 228 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE)); 229 } 230 #endif 231 232 if ((Ctx->dir0 == 0 && PetscAbsReal(t - 4.0) < 0.01) || (Ctx->dir0 == -1 && PetscAbsReal(t - 3.0) < 0.01)) { 233 SetVtols(Ctx->rank, Ctx->size, 1e-8, 1e-26, Ctx->vtol); // for better resolution of sin-event at t=5.0 234 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, Ctx->vtol)); 235 } 236 if (PetscAbsReal(t - 5.0) < 0.01) { 237 SetVtols(Ctx->rank, Ctx->size, 1e-8, 1e-8, Ctx->vtol); // back to normal 238 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, Ctx->vtol)); 239 } 240 241 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 242 PetscFunctionReturn(PETSC_SUCCESS); 243 } 244 245 // helper function to fill vtol[] 246 static inline void SetVtols(PetscMPIInt rank, PetscMPIInt size, PetscReal tol0, PetscReal tolsin, PetscReal *vtol) 247 { 248 PetscInt n = 0; 249 for (PetscInt i = -3; i <= 3; i++) 250 if (rank == (i + 3) % size) vtol[n++] = tol0; // pos-polynomials 251 for (PetscInt i = -3; i <= 3; i++) 252 if (rank == (i + 3) % size) vtol[n++] = tol0; // neg-polynomials 253 if (rank == 0) vtol[n++] = tolsin; // sin-event -- on rank-0 254 if (rank == size - 1) vtol[n++] = tol0; // cos-event -- on last rank 255 } 256 /*---------------------------------------------------------------------------------------------*/ 257 /* 258 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1, 259 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to 260 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed, 261 simply remove this option altogether. This will result in using the defaults 262 (which is PETSC_DECIDE). 263 */ 264 /*TEST 265 test: 266 suffix: pos1 267 output_file: output/ex5_pos1.out 268 args: -dir 1 -ts_event_dt_min 1e-6 269 args: -restart 1 270 args: -dtpost {{0 0.25}} 271 args: -ts_event_post_event_step 0.31 272 args: -ts_type rk 273 args: -ts_adapt_type {{none basic}} 274 nsize: 1 275 276 test: 277 suffix: pos4 278 output_file: output/ex5_pos4.out 279 args: -dir 1 -ts_event_dt_min 1e-6 -ts_dt 0.25 280 args: -restart 0 281 args: -dtpost 0 282 args: -ts_event_post_event_step -1 283 args: -ts_type {{beuler rk}} 284 args: -ts_adapt_type {{none basic}} 285 nsize: 4 286 filter: sort 287 filter_output: sort 288 289 test: 290 suffix: neu1 291 output_file: output/ex5_neu1.out 292 args: -dir 0 -ts_event_dt_min 1e-6 293 args: -restart 0 294 args: -dtpost {{0 0.25}} 295 args: -ts_event_post_event_step -1 296 args: -ts_type rk 297 args: -ts_adapt_type {{none basic}} 298 nsize: 1 299 300 test: 301 suffix: neu4 302 output_file: output/ex5_neu4.out 303 args: -dir 0 -ts_event_dt_min 1e-6 -ts_dt 0.25 304 args: -dtpost 0 305 args: -ts_event_post_event_step {{-1 0.29}} 306 args: -ts_event_post_event_second_step {{-1 0.31}} 307 args: -ts_type rk 308 args: -ts_adapt_type {{none basic}} 309 nsize: 4 310 filter: sort 311 filter_output: sort 312 313 test: 314 suffix: neg2 315 output_file: output/ex5_neg2.out 316 args: -dir -1 -ts_event_dt_min 1e-6 317 args: -restart 1 318 args: -dtpost {{0 0.25}} 319 args: -ts_event_post_event_step 0.31 320 args: -ts_type beuler 321 args: -ts_adapt_type {{none basic}} 322 nsize: 2 323 filter: sort 324 filter_output: sort 325 326 test: 327 suffix: neg4 328 output_file: output/ex5_neg4.out 329 args: -dir -1 -ts_event_dt_min 1e-6 -ts_dt 0.25 330 args: -restart 0 331 args: -dtpost 0 332 args: -ts_event_post_event_step -1 333 args: -ts_type {{beuler rk}} 334 args: -ts_adapt_type {{none basic}} 335 nsize: 4 336 filter: sort 337 filter_output: sort 338 339 TEST*/ 340