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 V*sin(pi*t), zeros = 1,...,10\n" 13 "After each event location the tolerance for the sin() event is multiplied by 4\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 "-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" 20 " if x == 0, nothing happens\n" 21 "-V {float}: scaling of the sin() event function; for small V this event is triggered by the function values, " 22 "for large V the event is triggered by the small step size\n" 23 "-change5 : flag to change the state vector at t=5 PostEvent\n"; 24 25 #define MAX_NFUNC 100 // max event functions per rank 26 #define MAX_NEV 5000 // max zero crossings for each rank 27 28 typedef struct { 29 PetscMPIInt rank, size; 30 PetscReal pi; 31 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 32 PetscReal evres[MAX_NEV]; // times of found zero-crossings 33 PetscInt evnum[MAX_NEV]; // number of zero-crossings at each time 34 PetscInt cnt; // counter 35 PetscBool flg; // flag for additional print in PostEvent 36 PetscBool restart; // flag for TSRestartStep() in PostEvent 37 PetscReal dtpost; // post-event step 38 PetscInt postcnt; // counter for PostEvent calls 39 PetscReal V; // vertical scaling for sin() 40 PetscReal vtol[MAX_NFUNC]; // vtol array, with extra storage 41 PetscBool change5; // flag to change the state vector at t=5 PostEvent 42 } AppCtx; 43 44 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 45 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 46 47 int main(int argc, char **argv) 48 { 49 TS ts; 50 Mat A; 51 Vec sol; 52 PetscInt n, dir0, m = 0; 53 PetscReal tol = 1e-7; 54 PetscInt dir[MAX_NFUNC], inds[2]; 55 PetscBool term[MAX_NFUNC]; 56 PetscScalar *x, vals[4]; 57 AppCtx ctx; 58 TSConvergedReason reason; 59 60 PetscFunctionBeginUser; 61 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 62 setbuf(stdout, NULL); 63 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 64 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 65 ctx.pi = PetscAcosReal(-1.0); 66 ctx.cnt = 0; 67 ctx.flg = PETSC_FALSE; 68 ctx.restart = PETSC_FALSE; 69 ctx.dtpost = 0; 70 ctx.postcnt = 0; 71 ctx.V = 1.0; 72 ctx.change5 = PETSC_FALSE; 73 74 // The linear problem has a 2*2 matrix. The matrix is constant 75 if (ctx.rank == 0) m = 2; 76 inds[0] = 0; 77 inds[1] = 1; 78 vals[0] = 0; 79 vals[1] = 0.2; 80 vals[2] = -0.2; 81 vals[3] = 0; 82 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A)); 83 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 84 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 85 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 86 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 87 88 PetscCall(MatCreateVecs(A, &sol, NULL)); 89 PetscCall(VecGetArray(sol, &x)); 90 if (ctx.rank == 0) { // initial conditions 91 x[0] = 0; // sin(0) 92 x[1] = 1; // cos(0) 93 } 94 PetscCall(VecRestoreArray(sol, &x)); 95 96 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 97 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 98 99 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 100 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL)); 101 102 PetscCall(TSSetTimeStep(ts, 0.1)); 103 PetscCall(TSSetType(ts, TSBEULER)); 104 PetscCall(TSSetMaxSteps(ts, 10000)); 105 PetscCall(TSSetMaxTime(ts, 10.0)); 106 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 107 PetscCall(TSSetFromOptions(ts)); 108 109 // Set the event handling 110 dir0 = 0; 111 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction 112 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 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(PetscOptionsGetReal(NULL, NULL, "-V", &ctx.V, NULL)); 116 PetscCall(PetscOptionsGetBool(NULL, NULL, "-change5", &ctx.change5, NULL)); // flag to change the state vector at t=5 PostEvent 117 118 n = 0; // event counter 119 if (ctx.rank == 0) { // first event -- on rank-0 120 ctx.vtol[n] = tol * 10; 121 dir[n] = dir0; 122 term[n++] = PETSC_FALSE; 123 } 124 if (ctx.rank == ctx.size - 1) { // second event (with termination) -- on last rank 125 ctx.vtol[n] = tol * 10; 126 dir[n] = dir0; 127 term[n++] = PETSC_TRUE; 128 } 129 if (ctx.rank == 1 % ctx.size) { // third event -- on rank = 1%ctx.size 130 ctx.vtol[n] = tol; 131 dir[n] = dir0; 132 term[n++] = PETSC_FALSE; 133 } 134 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 135 PetscCall(TSSetEventTolerances(ts, tol, ctx.vtol)); 136 137 // Solution 138 PetscCall(TSSolve(ts, sol)); 139 PetscCall(TSGetConvergedReason(ts, &reason)); 140 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "CONVERGED REASON: %" PetscInt_FMT " (TS_CONVERGED_EVENT == %" PetscInt_FMT ")\n", (PetscInt)reason, (PetscInt)TS_CONVERGED_EVENT)); 141 142 // The 3 columns printed are: [RANK] [num. of events at the given time] [time of event] 143 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])); 144 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 145 146 PetscCall(MatDestroy(&A)); 147 PetscCall(TSDestroy(&ts)); 148 PetscCall(VecDestroy(&sol)); 149 150 PetscCall(PetscFinalize()); 151 return 0; 152 } 153 154 /* 155 User callback for defining the event-functions 156 */ 157 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 158 { 159 PetscInt n = 0; 160 AppCtx *Ctx = (AppCtx *)ctx; 161 162 PetscFunctionBeginUser; 163 // for the test purposes, event-functions are defined based on t 164 // first event -- on rank-0 165 if (Ctx->rank == 0) { 166 if (t < 2.05) gval[n++] = 0.5 * (1 - PetscPowReal(t - 2.05, 12)); 167 else gval[n++] = 0.5; 168 } 169 170 // second event -- on last rank 171 if (Ctx->rank == Ctx->size - 1) { 172 if (t > 8.05) gval[n++] = 0.25 * (1 - PetscPowReal(t - 8.05, 12)); 173 else gval[n++] = 0.25; 174 } 175 176 // third event -- on rank = 1%ctx.size 177 if (Ctx->rank == 1 % Ctx->size) { gval[n++] = Ctx->V * PetscSinReal(Ctx->pi * t); } 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 /* 182 User callback for the post-event stuff 183 */ 184 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 185 { 186 PetscInt n = 0; 187 PetscScalar *x; 188 AppCtx *Ctx = (AppCtx *)ctx; 189 190 PetscFunctionBeginUser; 191 if (Ctx->flg) { 192 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 193 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 194 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 195 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 196 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 197 } 198 199 // change the state vector near t=5.0 200 if (PetscAbsReal(t - 5.0) < 0.01 && Ctx->change5) { 201 PetscCall(VecGetArray(U, &x)); 202 if (Ctx->rank == 0) x[1] = -x[1]; 203 PetscCall(VecRestoreArray(U, &x)); 204 } 205 206 // update vtol's 207 if (Ctx->rank == 0) n++; // first event -- on rank-0 208 if (Ctx->rank == Ctx->size - 1) n++; // second event -- on last rank 209 if (Ctx->rank == 1 % Ctx->size) { // third event -- on rank = 1%ctx.size 210 if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "vtol for sin: %g -> ", (double)Ctx->vtol[n])); 211 Ctx->vtol[n] *= 4; 212 if (PetscAbsReal(t - 5.0) < 0.01) Ctx->vtol[n] /= 100; // one-off decrease 213 if (Ctx->flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g\n", (double)Ctx->vtol[n])); 214 n++; 215 } 216 PetscCall(TSSetEventTolerances(ts, 0, Ctx->vtol)); 217 218 if (Ctx->cnt < MAX_NEV && nev_zero > 0) { 219 Ctx->evres[Ctx->cnt] = t; 220 Ctx->evnum[Ctx->cnt++] = nev_zero; 221 } 222 223 #ifdef NEW_VERSION 224 Ctx->postcnt++; // sync 225 if (Ctx->dtpost > 0) { 226 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 227 else PetscCall(TSSetPostEventStep(ts, 0)); 228 } 229 #endif 230 231 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 /*---------------------------------------------------------------------------------------------*/ 235 /*TEST 236 test: 237 suffix: V 238 output_file: output/ex3_V.out 239 args: -ts_type beuler 240 args: -ts_adapt_type basic 241 args: -V {{1e2 1e4 1e6 1e8}} 242 args: -ts_adapt_dt_min 1e-6 243 args: -change5 {{0 1}} 244 nsize: 1 245 246 test: 247 suffix: neu1 248 output_file: output/ex3_neu1.out 249 args: -dir 0 250 args: -V 1e5 251 args: -ts_adapt_dt_min 1e-6 252 args: -restart 1 253 args: -dtpost 0.24 254 args: -ts_event_post_event_step {{0 0.31}} 255 args: -ts_type {{beuler rk}} 256 args: -ts_adapt_type {{none basic}} 257 nsize: 1 258 259 test: 260 suffix: neu2 261 output_file: output/ex3_neu2.out 262 args: -dir 0 263 args: -V 1e5 264 args: -ts_adapt_dt_min 1e-6 265 args: -restart 1 266 args: -dtpost 0 267 args: -ts_event_post_event_step {{0 0.31}} 268 args: -ts_type {{beuler rk}} 269 args: -ts_adapt_type {{none basic}} 270 nsize: 2 271 filter: sort 272 filter_output: sort 273 274 test: 275 suffix: neu4 276 output_file: output/ex3_neu4.out 277 args: -dir 0 278 args: -V 1e5 279 args: -ts_adapt_dt_min 1e-6 280 args: -restart {{0 1}} 281 args: -dtpost 0.24 282 args: -ts_event_post_event_step 0.21 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: pos1 291 output_file: output/ex3_pos1.out 292 args: -dir 1 293 args: -V 1e5 294 args: -ts_adapt_dt_min 1e-6 295 args: -restart 0 296 args: -dtpost {{0 0.24}} 297 args: -ts_event_post_event_step 0 298 args: -ts_type {{beuler rk}} 299 args: -ts_adapt_type {{none basic}} 300 nsize: 1 301 302 test: 303 suffix: pos2 304 output_file: output/ex3_pos2.out 305 args: -dir 1 306 args: -V 1e5 307 args: -ts_adapt_dt_min 1e-6 308 args: -restart 1 309 args: -dtpost {{0 0.24}} 310 args: -ts_event_post_event_step 0 311 args: -ts_type {{beuler rk}} 312 args: -ts_adapt_type {{none basic}} 313 nsize: 2 314 filter: sort 315 filter_output: sort 316 317 test: 318 suffix: pos4 319 output_file: output/ex3_pos4.out 320 args: -dir 1 321 args: -V 1e5 322 args: -ts_adapt_dt_min 1e-6 323 args: -restart 0 324 args: -dtpost 0 325 args: -ts_event_post_event_step {{0 0.32}} 326 args: -ts_type {{beuler rk}} 327 args: -ts_adapt_type {{none basic}} 328 nsize: 4 329 filter: sort 330 filter_output: sort 331 332 test: 333 suffix: neg1 334 output_file: output/ex3_neg1.out 335 args: -dir -1 336 args: -V 1e5 337 args: -ts_adapt_dt_min 1e-6 338 args: -restart 1 339 args: -dtpost {{0 0.24}} 340 args: -ts_event_post_event_step 0 341 args: -ts_type {{beuler rk}} 342 args: -ts_adapt_type {{none basic}} 343 nsize: 1 344 345 test: 346 suffix: neg2 347 output_file: output/ex3_neg2.out 348 args: -dir -1 349 args: -V 1e5 350 args: -ts_adapt_dt_min 1e-6 351 args: -restart 0 352 args: -dtpost {{0 0.24}} 353 args: -ts_event_post_event_step 0 354 args: -ts_type {{beuler rk}} 355 args: -ts_adapt_type {{none basic}} 356 nsize: 2 357 filter: sort 358 filter_output: sort 359 360 test: 361 suffix: neg4 362 output_file: output/ex3_neg4.out 363 args: -dir -1 364 args: -V 1e5 365 args: -ts_adapt_dt_min 1e-6 366 args: -restart 0 367 args: -dtpost {{0 0.24}} 368 args: -ts_event_post_event_step 0.3 369 args: -ts_type {{beuler rk}} 370 args: -ts_adapt_type {{none basic}} 371 nsize: 4 372 filter: sort 373 filter_output: sort 374 TEST*/ 375