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 event function = sin(pi*t), zeros = 1,...,10\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, (char *)0, 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 113 for (PetscInt i = 1; i < MAX_NEV; i++) { 114 if (i % 2 == 1 && dir0 <= 0) ctx.ref[ctx.cntref++] = i; 115 if (i % 2 == 0 && dir0 >= 0) ctx.ref[ctx.cntref++] = i; 116 } 117 } 118 if (ctx.cntref > 0) PetscCall(PetscSortReal(ctx.cntref, ctx.ref)); 119 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 120 121 // Solution 122 PetscCall(TSSolve(ts, sol)); 123 124 // The 4 columns printed are: [RANK] [time of event] [error w.r.t. reference] ["pass"/"fail"] 125 for (PetscInt j = 0; j < ctx.cnt; j++) { 126 PetscReal err = 10.0; 127 if (j < ctx.cntref) err = PetscAbsReal(ctx.evres[j] - ctx.ref[j]); 128 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")); 129 } 130 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 131 132 PetscCall(MatDestroy(&A)); 133 PetscCall(TSDestroy(&ts)); 134 PetscCall(VecDestroy(&sol)); 135 136 PetscCall(PetscFinalize()); 137 return 0; 138 } 139 140 /* 141 User callback for defining the event-functions 142 */ 143 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 144 { 145 PetscInt n = 0; 146 AppCtx *Ctx = (AppCtx *)ctx; 147 148 PetscFunctionBeginUser; 149 // for the test purposes, event-functions are defined based on t 150 // first event -- on rank-0 151 if (Ctx->rank == 0) { gval[n++] = PetscSinReal(Ctx->pi * t); } 152 PetscFunctionReturn(PETSC_SUCCESS); 153 } 154 155 /* 156 User callback for the post-event stuff 157 */ 158 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 159 { 160 AppCtx *Ctx = (AppCtx *)ctx; 161 162 PetscFunctionBeginUser; 163 if (Ctx->flg) { 164 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 165 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 166 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 167 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 168 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 169 } 170 171 if (Ctx->cnt + nev_zero < MAX_NEV) 172 for (PetscInt i = 0; i < nev_zero; i++) Ctx->evres[Ctx->cnt++] = t; // save the repeating zeros separately for easier/unified testing 173 174 #ifdef NEW_VERSION 175 Ctx->postcnt++; // sync 176 if (Ctx->dtpost > 0) { 177 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 178 else PetscCall(TSSetPostEventStep(ts, PETSC_DECIDE)); 179 } 180 #endif 181 182 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 183 PetscFunctionReturn(PETSC_SUCCESS); 184 } 185 /*---------------------------------------------------------------------------------------------*/ 186 /* 187 Note, in the tests below, -ts_event_post_event_step is occasionally set to -1, 188 which corresponds to PETSC_DECIDE in the API. It is not a very good practice to 189 explicitly specify -1 in this option. Rather, if PETSC_DECIDE behaviour is needed, 190 simply remove this option altogether. This will result in using the defaults 191 (which is PETSC_DECIDE). 192 */ 193 /*TEST 194 test: 195 suffix: 0 196 requires: !single 197 output_file: output/ex1sin_0.out 198 args: -dir 0 199 args: -restart 0 200 args: -dtpost {{0 0.25}} 201 args: -ts_event_post_event_step -1 202 args: -ts_type {{beuler rk}} 203 args: -ts_adapt_type {{none basic}} 204 nsize: 1 205 206 test: 207 suffix: p 208 requires: !single 209 output_file: output/ex1sin_p.out 210 args: -dir 1 211 args: -restart 1 212 args: -dtpost {{0 0.31}} 213 args: -ts_event_post_event_step {{-1 0.25}} 214 args: -ts_type rk 215 args: -ts_adapt_type {{none basic}} 216 nsize: 2 217 filter: sort 218 filter_output: sort 219 220 test: 221 suffix: n 222 requires: !single 223 output_file: output/ex1sin_n.out 224 args: -dir -1 225 args: -restart {{0 1}} 226 args: -dtpost {{0 0.25}} 227 args: -ts_event_post_event_step 0.31 228 args: -ts_type {{beuler rk}} 229 args: -ts_adapt_type basic 230 nsize: 4 231 filter: sort 232 filter_output: sort 233 234 test: 235 suffix: 0single 236 requires: single 237 output_file: output/ex1sin_0.out 238 args: -dir 0 -ts_event_dt_min 1e-6 -errtol 5e-5 239 args: -restart 1 240 args: -dtpost {{0 0.25}} 241 args: -ts_event_post_event_step {{-1 0.31}} 242 args: -ts_type beuler 243 args: -ts_adapt_type {{none basic}} 244 nsize: 4 245 filter: sort 246 filter_output: sort 247 248 test: 249 suffix: psingle 250 requires: single 251 output_file: output/ex1sin_p.out 252 args: -dir 1 -ts_event_dt_min 1e-6 -errtol 5e-5 253 args: -restart 0 254 args: -dtpost {{0 0.31}} 255 args: -ts_event_post_event_step 0.25 256 args: -ts_type {{beuler rk}} 257 args: -ts_adapt_type {{none basic}} 258 nsize: 1 259 260 test: 261 suffix: nsingle 262 requires: single 263 output_file: output/ex1sin_n.out 264 args: -dir -1 -ts_event_dt_min 1e-6 -errtol 5e-5 265 args: -restart 1 266 args: -dtpost {{0 0.25}} 267 args: -ts_event_post_event_step -1 268 args: -ts_type {{beuler rk}} 269 args: -ts_adapt_type {{none basic}} 270 nsize: 2 271 filter: sort 272 filter_output: sort 273 TEST*/ 274