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 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 "-restart : flag for TSRestartStep() in PostEvent\n" 14 "-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" 15 " if x == 0, nothing happens\n"; 16 17 #define MAX_NFUNC 100 // max event functions per rank 18 #define MAX_NEV 5000 // max zero crossings for each rank 19 20 typedef struct { 21 PetscMPIInt rank, size; 22 PetscReal pi; 23 PetscReal fvals[MAX_NFUNC]; // helper array for reporting the residuals 24 PetscReal evres[MAX_NEV]; // times of found zero-crossings 25 PetscInt evnum[MAX_NEV]; // number of zero-crossings at each time 26 PetscInt cnt; // counter 27 PetscBool flg; // flag for additional print in PostEvent 28 PetscBool restart; // flag for TSRestartStep() in PostEvent 29 PetscReal dtpost; // post-event step 30 PetscInt postcnt; // counter for PostEvent calls 31 } AppCtx; 32 33 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx); 34 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx); 35 36 int main(int argc, char **argv) 37 { 38 TS ts; 39 Mat A; 40 Vec sol; 41 PetscInt n, dir0, m = 0; 42 PetscInt dir[MAX_NFUNC], inds[2]; 43 PetscBool term[MAX_NFUNC]; 44 PetscScalar *x, vals[4]; 45 AppCtx ctx; 46 47 PetscFunctionBeginUser; 48 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 49 setbuf(stdout, NULL); 50 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank)); 51 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size)); 52 ctx.pi = PetscAcosReal(-1.0); 53 ctx.cnt = 0; 54 ctx.flg = PETSC_FALSE; 55 ctx.restart = PETSC_FALSE; 56 ctx.dtpost = 0; 57 ctx.postcnt = 0; 58 59 // The linear problem has a 2*2 matrix. The matrix is constant 60 if (ctx.rank == 0) m = 2; 61 inds[0] = 0; 62 inds[1] = 1; 63 vals[0] = 0; 64 vals[1] = 0.2; 65 vals[2] = -0.2; 66 vals[3] = 0; 67 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, m, m, PETSC_DETERMINE, PETSC_DETERMINE, 2, NULL, 0, NULL, &A)); 68 PetscCall(MatSetValues(A, m, inds, m, inds, vals, INSERT_VALUES)); 69 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 70 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 71 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 72 73 PetscCall(MatCreateVecs(A, &sol, NULL)); 74 PetscCall(VecGetArray(sol, &x)); 75 if (ctx.rank == 0) { // initial conditions 76 x[0] = 0; // sin(0) 77 x[1] = 1; // cos(0) 78 } 79 PetscCall(VecRestoreArray(sol, &x)); 80 81 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 82 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 83 84 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, NULL)); 85 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, NULL)); 86 87 PetscCall(TSSetTimeStep(ts, 0.1)); 88 PetscCall(TSSetType(ts, TSBEULER)); 89 PetscCall(TSSetMaxSteps(ts, 10000)); 90 PetscCall(TSSetMaxTime(ts, 10.0)); 91 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 92 PetscCall(TSSetFromOptions(ts)); 93 94 // Set the event handling 95 dir0 = 0; 96 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &dir0, NULL)); // desired zero-crossing direction 97 PetscCall(PetscOptionsHasName(NULL, NULL, "-flg", &ctx.flg)); // flag for additional output 98 PetscCall(PetscOptionsGetBool(NULL, NULL, "-restart", &ctx.restart, NULL)); // flag for TSRestartStep() 99 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dtpost", &ctx.dtpost, NULL)); // post-event step 100 101 n = 0; // event counter 102 if (ctx.rank == 0) { // first event -- on rank-0 103 dir[n] = dir0; 104 term[n++] = PETSC_FALSE; 105 } 106 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 107 108 // Solution 109 PetscCall(TSSolve(ts, sol)); 110 111 // The 3 columns printed are: [RANK] [num. of events at the given time] [time of event] 112 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])); 113 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 114 115 PetscCall(MatDestroy(&A)); 116 PetscCall(TSDestroy(&ts)); 117 PetscCall(VecDestroy(&sol)); 118 119 PetscCall(PetscFinalize()); 120 return 0; 121 } 122 123 /* 124 User callback for defining the event-functions 125 */ 126 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 127 { 128 PetscInt n = 0; 129 AppCtx *Ctx = (AppCtx *)ctx; 130 131 PetscFunctionBeginUser; 132 // for the test purposes, event-functions are defined based on t 133 // first event -- on rank-0 134 if (Ctx->rank == 0) { gval[n++] = PetscSinReal(Ctx->pi * t); } 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 /* 139 User callback for the post-event stuff 140 */ 141 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 142 { 143 AppCtx *Ctx = (AppCtx *)ctx; 144 145 PetscFunctionBeginUser; 146 if (Ctx->flg) { 147 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 148 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 149 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 150 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 151 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 152 } 153 154 if (Ctx->cnt < MAX_NEV && nev_zero > 0) { 155 Ctx->evres[Ctx->cnt] = t; 156 Ctx->evnum[Ctx->cnt++] = nev_zero; 157 } 158 159 #ifdef NEW_VERSION 160 Ctx->postcnt++; // sync 161 if (Ctx->dtpost > 0) { 162 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 163 else PetscCall(TSSetPostEventStep(ts, 0)); 164 } 165 #endif 166 167 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170 /*---------------------------------------------------------------------------------------------*/ 171 /*TEST 172 test: 173 suffix: 0 174 output_file: output/ex1sin_0.out 175 args: -dir 0 176 args: -restart {{0 1}} 177 args: -dtpost {{0 0.25}} 178 args: -ts_event_post_event_step {{0 0.31}} 179 args: -ts_type {{beuler rk}} 180 args: -ts_adapt_type {{none basic}} 181 nsize: {{1 4}} 182 filter: sort 183 filter_output: sort 184 185 test: 186 suffix: p 187 output_file: output/ex1sin_p.out 188 args: -dir 1 189 args: -restart {{0 1}} 190 args: -dtpost {{0 0.31}} 191 args: -ts_event_post_event_step {{0 0.25}} 192 args: -ts_type {{beuler rk}} 193 args: -ts_adapt_type {{none basic}} 194 nsize: {{1 2}} 195 filter: sort 196 filter_output: sort 197 198 test: 199 suffix: n 200 output_file: output/ex1sin_n.out 201 args: -dir -1 202 args: -restart {{0 1}} 203 args: -dtpost {{0 0.25}} 204 args: -ts_event_post_event_step {{0 0.31}} 205 args: -ts_type {{beuler rk}} 206 args: -ts_adapt_type {{none basic}} 207 nsize: {{1 4}} 208 filter: sort 209 filter_output: sort 210 TEST*/ 211