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 or several event functions (on rank-0)\n" 10 "This program is mostly intended to test the Anderson-Bjorck iteration\n" 11 "Options:\n" 12 "-dir d : zero-crossing direction for events\n" 13 "-flg : additional output in Postevent\n" 14 "-restart : flag for TSRestartStep() in PostEvent\n" 15 "-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" 16 " if x == 0, nothing happens\n" 17 "-func F : selects the event function [0, ..., 11], if F == -1 (default) is set, all event functions are taken\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 PetscInt evnum[MAX_NEV]; // number of zero-crossings at each time 28 PetscInt cnt; // counter 29 PetscBool flg; // flag for additional print in PostEvent 30 PetscBool restart; // flag for TSRestartStep() in PostEvent 31 PetscReal dtpost; // post-event step 32 PetscInt postcnt; // counter for PostEvent calls 33 PetscInt F; // event-function index 34 PetscInt Fnum; // total available event functions 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.flg = PETSC_FALSE; 59 ctx.restart = PETSC_FALSE; 60 ctx.dtpost = 0; 61 ctx.postcnt = 0; 62 ctx.F = -1; 63 ctx.Fnum = 12; 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(TSSetTime(ts, 0.03)); 94 PetscCall(TSSetTimeStep(ts, 0.1)); 95 PetscCall(TSSetType(ts, TSBEULER)); 96 PetscCall(TSSetMaxSteps(ts, 10000)); 97 PetscCall(TSSetMaxTime(ts, 4.0)); 98 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 99 PetscCall(TSSetFromOptions(ts)); 100 101 // Set the event handling 102 dir0 = 0; 103 PetscCall(PetscOptionsGetInt(NULL, NULL, "-dir", &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 PetscCall(PetscOptionsGetInt(NULL, NULL, "-F", &ctx.F, NULL)); // event-function index 108 PetscCheck(ctx.F >= -1 && ctx.F < ctx.Fnum, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_OUTOFRANGE, "Value of 'F' is out of range"); 109 110 n = 0; // event counter 111 if (ctx.rank == 0) { // all events -- on rank-0 112 if (ctx.F == -1) 113 for (n = 0; n < ctx.Fnum; n++) { // all event-functions 114 dir[n] = dir0; 115 term[n] = PETSC_FALSE; 116 } 117 else { // single event-function 118 dir[n] = dir0; 119 term[n++] = PETSC_FALSE; 120 } 121 } 122 PetscCall(TSSetEventHandler(ts, n, dir, term, EventFunction, Postevent, &ctx)); 123 124 // Solution 125 PetscCall(TSSolve(ts, sol)); 126 127 // The 3 columns printed are: [RANK] [num. of events at the given time] [time of event] 128 for (PetscInt j = 0; j < ctx.cnt; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%d\t%" PetscInt_FMT "\t%.5g\n", ctx.rank, ctx.evnum[j], (double)ctx.evres[j])); 129 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 130 131 PetscCall(MatDestroy(&A)); 132 PetscCall(TSDestroy(&ts)); 133 PetscCall(VecDestroy(&sol)); 134 135 PetscCall(PetscFinalize()); 136 return 0; 137 } 138 139 /* 140 User callback for defining the event-functions 141 */ 142 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal gval[], void *ctx) 143 { 144 PetscInt n = 0; 145 AppCtx *Ctx = (AppCtx *)ctx; 146 147 PetscFunctionBeginUser; 148 // for the test purposes, event-functions are defined based on t 149 // all events -- on rank-0 150 if (Ctx->rank == 0) { 151 if (Ctx->F == 0 || Ctx->F == -1) gval[n++] = PetscSinReal(Ctx->pi * t) / Ctx->pi; // FUNC-0, roots 1, 2, 3, 4 152 if (Ctx->F == 1 || Ctx->F == -1) gval[n++] = PetscLogReal(t); // FUNC-2, root 1 153 if (Ctx->F == 2 || Ctx->F == -1) { // FUNC-3, root 1 154 if (t < 2) gval[n++] = (1 - PetscPowReal(t - 2, 12)) / 12.0; 155 else gval[n++] = 1 / 12.0; 156 } 157 if (Ctx->F == 3 || Ctx->F == -1) gval[n++] = t - PetscExpReal(PetscSinReal(t)) + 1; // FUNC-5, root 1.69681 158 if (Ctx->F == 4 || Ctx->F == -1) gval[n++] = (1e10 * PetscPowReal(t, 1 / t) - 1) / 100; // FUNC-6, root 0.1 159 if (Ctx->F == 5 || Ctx->F == -1) gval[n++] = PetscLogReal(t - 0.02) * PetscLogReal(t - 0.02) * PetscSignReal(t - 1.02) * 1e7; // FUNC-7, root 1.02 160 if (Ctx->F == 6 || Ctx->F == -1) gval[n++] = 4 * PetscCosReal(t) - PetscExpReal(t); // FUNC-14, root 0.904788 161 if (Ctx->F == 7 || Ctx->F == -1) gval[n++] = (20.0 * t - 1) / (19.0 * t) / 10; // FUNC-15, root 0.05 162 if (Ctx->F == 8 || Ctx->F == -1) gval[n++] = ((t - 1) * PetscExpReal(-20 * t) + PetscPowReal(t, 20)) * 1e4; // FUNC-16, root 0.552 163 if (Ctx->F == 9 || Ctx->F == -1) gval[n++] = (t * t * (t * t / 3.0 + PetscSqrtReal(2.0) * PetscSinReal(t)) - PetscSqrtReal(3.0) / 18) * 10; // FUNC-17, root 0.399 164 if (Ctx->F == 10 || Ctx->F == -1) gval[n++] = ((t * t + 1) * PetscSinReal(t) - PetscExpReal(PetscSqrtReal(t)) * (t - 1) * (t * t - 5)) / 10; // FUNC-18, roots 0.87, 2.388 165 if (Ctx->F == 11 || Ctx->F == -1) gval[n++] = 2 * t - 5; // FUNC-21, root 2.5 166 } 167 PetscFunctionReturn(PETSC_SUCCESS); 168 } 169 170 /* 171 User callback for the post-event stuff 172 */ 173 PetscErrorCode Postevent(TS ts, PetscInt nev_zero, PetscInt evs_zero[], PetscReal t, Vec U, PetscBool fwd, void *ctx) 174 { 175 AppCtx *Ctx = (AppCtx *)ctx; 176 177 PetscFunctionBeginUser; 178 if (Ctx->flg) { 179 PetscCallBack("EventFunction", EventFunction(ts, t, U, Ctx->fvals, ctx)); 180 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] At t = %20.16g : %" PetscInt_FMT " events triggered, fvalues =", Ctx->rank, (double)t, nev_zero)); 181 for (PetscInt j = 0; j < nev_zero; j++) PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\t%g", (double)Ctx->fvals[evs_zero[j]])); 182 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n")); 183 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 184 } 185 186 if (Ctx->cnt + nev_zero < MAX_NEV) { 187 for (PetscInt i = 0; i < nev_zero; i++) { // save the repeating zeros separately for easier/unified testing 188 Ctx->evres[Ctx->cnt] = t; 189 Ctx->evnum[Ctx->cnt++] = 1; 190 } 191 } 192 193 #ifdef NEW_VERSION 194 Ctx->postcnt++; // sync 195 if (Ctx->dtpost > 0) { 196 if (Ctx->postcnt % 2 == 0) PetscCall(TSSetPostEventStep(ts, Ctx->dtpost)); 197 else PetscCall(TSSetPostEventStep(ts, 0)); 198 } 199 #endif 200 201 if (Ctx->restart) PetscCall(TSRestartStep(ts)); 202 PetscFunctionReturn(PETSC_SUCCESS); 203 } 204 /*---------------------------------------------------------------------------------------------*/ 205 /*TEST 206 test: 207 suffix: 0 208 output_file: output/ex4_0.out 209 args: -dir 0 210 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6 211 args: -ts_dt 0.4 212 args: -restart 1 213 args: -ts_event_tol {{1e-8 1e-15}} 214 args: -ts_adapt_type {{none basic}} 215 args: -dtpost {{0 0.25}} 216 args: -ts_event_post_event_step {{0 0.35}} 217 args: -ts_type {{beuler rk}} 218 nsize: {{1 4}} 219 filter: sort 220 filter_output: sort 221 222 test: 223 suffix: F7 224 output_file: output/ex4_F7.out 225 args: -dir 0 226 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6 227 args: -ts_dt 0.4 228 args: -F 7 229 args: -ts_event_tol {{1e-8 1e-15}} 230 args: -ts_adapt_type {{none basic}} 231 args: -ts_type {{beuler rk}} 232 nsize: 1 233 234 test: 235 suffix: F7revisit 236 output_file: output/ex4_F7revisit.out 237 args: -ts_event_monitor -F 7 -ts_dt 0.04 -ts_event_dt_min 0.016 238 nsize: 1 239 240 test: 241 suffix: pos 242 output_file: output/ex4_pos.out 243 args: -dir 1 244 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6 245 args: -ts_dt 0.4 246 args: -restart 0 247 args: -ts_event_tol {{1e-8 1e-15}} 248 args: -ts_adapt_type {{none basic}} 249 args: -dtpost {{0 0.25}} 250 args: -ts_event_post_event_step {{0 0.35}} 251 args: -ts_type {{beuler rk}} 252 nsize: {{1 4}} 253 filter: sort 254 filter_output: sort 255 256 test: 257 suffix: neg 258 output_file: output/ex4_neg.out 259 args: -dir -1 260 args: -ts_adapt_dt_min 1e-10 -ts_event_dt_min 1e-6 261 args: -ts_dt 0.4 262 args: -restart 1 263 args: -ts_event_tol {{1e-8 1e-15}} 264 args: -ts_adapt_type {{none basic}} 265 args: -dtpost {{0 0.25}} 266 args: -ts_event_post_event_step {{0 0.35}} 267 args: -ts_type {{beuler rk}} 268 nsize: {{1 4}} 269 filter: sort 270 filter_output: sort 271 TEST*/ 272