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 "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 "-restart : flag for TSRestartStep() in PostEvent\n" 18 "-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" 19 " if x == 0, nothing happens\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 PetscInt evnum[MAX_NEV]; // number of zero-crossings at each time 30 PetscInt cnt; // counter 31 PetscBool flg; // flag for additional print in PostEvent 32 PetscBool restart; // flag for TSRestartStep() in PostEvent 33 PetscReal dtpost; // post-event step 34 PetscInt postcnt; // counter for PostEvent calls 35 PetscReal vtol[MAX_NFUNC]; // vtol array, with extra storage 36 PetscInt dir0; // desired zero-crossing direction 37 } AppCtx; 38 39 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 40 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 41 static inline void SetVtols(PetscMPIInt rank, PetscMPIInt size, PetscReal tol0, PetscReal tolsin, PetscReal *vtol); // helper function to fill vtol[] 42 43 int main(int argc, char **argv) 44 { 45 TS ts; 46 Mat A; 47 Vec sol; 48 PetscInt n, m = 0; 49 PetscInt dir[MAX_NFUNC], inds[2]; 50 PetscBool term[MAX_NFUNC]; 51 PetscScalar *x, vals[4]; 52 AppCtx ctx; 53 54 PetscFunctionBeginUser; 55 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 56 setbuf(stdout, NULL); 57 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 58 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 59 ctx.pi = PetscAcosReal(-1.0); 60 ctx.cnt = 0; 61 ctx.flg = PETSC_FALSE; 62 ctx.restart = PETSC_FALSE; 63 ctx.dtpost = 0; 64 ctx.postcnt = 0; 65 66 // The linear problem has a 2*2 matrix. The matrix is constant 67 if (ctx.rank == 0) m = 2; 68 inds[0] = 0; 69 inds[1] = 1; 70 vals[0] = 0; 71 vals[1] = 0.2; 72 vals[2] = -0.2; 73 vals[3] = 0; 74 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A)); 75 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 76 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 77 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 78 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 79 80 PetscCall(MatCreateVecs(A, &sol, NULL)); 81 PetscCall(VecGetArray(sol, &x)); 82 if (ctx.rank == 0) { // initial conditions 83 x[0] = 0; // sin(0) 84 x[1] = 1; // cos(0) 85 } 86 PetscCall(VecRestoreArray(sol, &x)); 87 88 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 89 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 90 91 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 92 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL)); 93 94 PetscCall(TSSetTimeStep(ts, 0.1)); 95 PetscCall(TSSetType(ts, TSBEULER)); 96 PetscCall(TSSetMaxSteps(ts, 10000)); 97 PetscCall(TSSetMaxTime(ts, 10.0)); 98 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 99 PetscCall(TSSetFromOptions(ts)); 100 101 // Set the event handling 102 ctx.dir0 = 0; 103 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &ctx.dir0, NULL)); // desired zero-crossing direction 104 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 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 for (PetscInt i = -3; i <= 3; i++) { // pos-polynomials 110 if (ctx.rank == (i + 3) % ctx.size) { 111 dir[n] = ctx.dir0; 112 term[n++] = PETSC_FALSE; 113 } 114 } 115 for (PetscInt i = -3; i <= 3; i++) { // neg-polynomials 116 if (ctx.rank == (i + 3) % ctx.size) { 117 dir[n] = ctx.dir0; 118 term[n++] = PETSC_FALSE; 119 } 120 } 121 if (ctx.rank == 0) { // sin-event -- on rank-0 122 dir[n] = ctx.dir0; 123 term[n++] = PETSC_FALSE; 124 } 125 if (ctx.rank == ctx.size - 1) { // cos-event -- on last rank 126 dir[n] = ctx.dir0; 127 term[n++] = PETSC_FALSE; 128 } 129 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 130 SetVtols(ctx.rank, ctx.size, 1e-8, 1e-8, ctx.vtol); 131 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, ctx.vtol)); 132 133 // Solution 134 PetscCall(TSSolve(ts, sol)); 135 136 // The 3 columns printed are: [RANK] [num. of events at the given time] [time of event] 137 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])); 138 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 139 140 PetscCall(MatDestroy(&A)); 141 PetscCall(TSDestroy(&ts)); 142 PetscCall(VecDestroy(&sol)); 143 144 PetscCall(PetscFinalize()); 145 return 0; 146 } 147 148 /* 149 User callback for defining the event-functions 150 */ 151 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 152 { 153 PetscInt n = 0; 154 AppCtx *Ctx = (AppCtx *)ctx; 155 PetscReal P; 156 157 PetscFunctionBeginUser; 158 // for the test purposes, event-functions are defined based on t 159 for (PetscInt i = -3; i <= 3; i++) { // pos-polynomials 160 if (Ctx->rank == (i + 3) % Ctx->size) { 161 P = PetscPowReal(2.0, i); 162 if (t < 2 + P) gval[n++] = 1 - PetscPowReal(2 + P - t, i + 5); 163 else gval[n++] = 1; 164 } 165 } 166 for (PetscInt i = -3; i <= 3; i++) { // neg-polynomials 167 if (Ctx->rank == (i + 3) % Ctx->size) { 168 P = PetscPowReal(2.0, i); 169 if (t > 8 - P) gval[n++] = 1 - PetscPowReal(t - 8 + P, i + 5); 170 else gval[n++] = 1; 171 } 172 } 173 if (Ctx->rank == 0) gval[n++] = (t - 5) * (t - 5) * PetscSinReal(Ctx->pi * t); // sin-event -- on rank-0 174 if (Ctx->rank == Ctx->size - 1) gval[n++] = 0.5 * PetscCosReal(Ctx->pi * t); // cos-event -- on last rank 175 PetscFunctionReturn(PETSC_SUCCESS); 176 } 177 178 /* 179 User callback for the post-event stuff 180 */ 181 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 182 { 183 AppCtx *Ctx = (AppCtx *)ctx; 184 185 PetscFunctionBeginUser; 186 if (Ctx->flg) { 187 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 188 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 189 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 190 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 191 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 192 } 193 194 if (Ctx->cnt + nev_zero < MAX_NEV) { 195 for (PetscInt i = 0; i < nev_zero; i++) { // save the repeating zeros separately for easier/unified testing 196 Ctx->evres[Ctx->cnt] = t; 197 Ctx->evnum[Ctx->cnt++] = 1; 198 } 199 } 200 201 #ifdef NEW_VERSION 202 Ctx->postcnt++; // sync 203 if (Ctx->dtpost > 0) { 204 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 205 else PetscCall(TSSetPostEventStep(ts, 0)); 206 } 207 #endif 208 209 if ((Ctx->dir0 == 0 && PetscAbsReal(t - 4.0) < 0.01) || (Ctx->dir0 == -1 && PetscAbsReal(t - 3.0) < 0.01)) { 210 SetVtols(Ctx->rank, Ctx->size, 1e-8, 1e-26, Ctx->vtol); // for better resolution of sin-event at t=5.0 211 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, Ctx->vtol)); 212 } 213 if (PetscAbsReal(t - 5.0) < 0.01) { 214 SetVtols(Ctx->rank, Ctx->size, 1e-8, 1e-8, Ctx->vtol); // back to normal 215 PetscCall(TSSetEventTolerances(ts, PETSC_DECIDE, Ctx->vtol)); 216 } 217 218 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 219 PetscFunctionReturn(PETSC_SUCCESS); 220 } 221 222 // helper function to fill vtol[] 223 static inline void SetVtols(PetscMPIInt rank, PetscMPIInt size, PetscReal tol0, PetscReal tolsin, PetscReal *vtol) 224 { 225 PetscInt n = 0; 226 for (PetscInt i = -3; i <= 3; i++) 227 if (rank == (i + 3) % size) vtol[n++] = tol0; // pos-polynomials 228 for (PetscInt i = -3; i <= 3; i++) 229 if (rank == (i + 3) % size) vtol[n++] = tol0; // neg-polynomials 230 if (rank == 0) vtol[n++] = tolsin; // sin-event -- on rank-0 231 if (rank == size - 1) vtol[n++] = tol0; // cos-event -- on last rank 232 } 233 /*---------------------------------------------------------------------------------------------*/ 234 /*TEST 235 236 test: 237 suffix: pos1 238 output_file: output/ex5_pos1.out 239 args: -dir 1 -ts_event_dt_min 1e-6 240 args: -restart {{0 1}} 241 args: -dtpost {{0 0.25}} 242 args: -ts_event_post_event_step {{0 0.31}} 243 args: -ts_type {{beuler rk}} 244 args: -ts_adapt_type {{none basic}} 245 nsize: 1 246 247 test: 248 suffix: pos4 249 output_file: output/ex5_pos4.out 250 args: -dir 1 -ts_event_dt_min 1e-6 251 args: -restart {{0 1}} 252 args: -dtpost {{0 0.25}} 253 args: -ts_event_post_event_step {{0 0.31}} 254 args: -ts_type {{beuler rk}} 255 args: -ts_adapt_type {{none basic}} 256 nsize: 4 257 filter: sort 258 filter_output: sort 259 260 test: 261 suffix: neu1 262 output_file: output/ex5_neu1.out 263 args: -dir 0 -ts_event_dt_min 1e-6 264 args: -restart {{0 1}} 265 args: -dtpost {{0 0.25}} 266 args: -ts_event_post_event_step {{0 0.31}} 267 args: -ts_type {{beuler rk}} 268 args: -ts_adapt_type {{none basic}} 269 nsize: 1 270 271 test: 272 suffix: neu4 273 output_file: output/ex5_neu4.out 274 args: -dir 0 -ts_event_dt_min 1e-6 275 args: -restart {{0 1}} 276 args: -dtpost {{0 0.25}} 277 args: -ts_event_post_event_step {{0 0.31}} 278 args: -ts_type {{beuler rk}} 279 args: -ts_adapt_type {{none basic}} 280 nsize: 4 281 filter: sort 282 filter_output: sort 283 284 test: 285 suffix: neg2 286 output_file: output/ex5_neg2.out 287 args: -dir -1 -ts_event_dt_min 1e-6 288 args: -restart {{0 1}} 289 args: -dtpost {{0 0.25}} 290 args: -ts_event_post_event_step {{0 0.31}} 291 args: -ts_type {{beuler rk}} 292 args: -ts_adapt_type {{none basic}} 293 nsize: 2 294 filter: sort 295 filter_output: sort 296 297 test: 298 suffix: neg4 299 output_file: output/ex5_neg4.out 300 args: -dir -1 -ts_event_dt_min 1e-6 301 args: -restart {{0 1}} 302 args: -dtpost {{0 0.25}} 303 args: -ts_event_post_event_step {{0 0.31}} 304 args: -ts_type {{beuler rk}} 305 args: -ts_adapt_type {{none basic}} 306 nsize: 4 307 filter: sort 308 filter_output: sort 309 310 TEST*/ 311