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