1 static char help[] = "Solves the trivial ODE du/dt = 1, u(0) = 0. \n\n"; 2 3 #include <petscts.h> 4 #include <petscpc.h> 5 6 static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 7 static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *); 8 9 static PetscErrorCode PreStep(TS); 10 static PetscErrorCode PostStep(TS); 11 static PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *); 12 static PetscErrorCode Event(TS, PetscReal, Vec, PetscReal *, void *); 13 static PetscErrorCode PostEvent(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *); 14 static PetscErrorCode TransferSetUp(TS, PetscInt, PetscReal, Vec, PetscBool *, void *); 15 static PetscErrorCode Transfer(TS, PetscInt, Vec[], Vec[], void *); 16 17 int main(int argc, char **argv) 18 { 19 TS ts; 20 PetscInt n, ntransfer[] = {2, 2}; 21 const PetscInt n_end = 11; 22 PetscReal t; 23 const PetscReal t_end = 11; 24 Vec x; 25 Vec f; 26 Mat A; 27 28 PetscFunctionBeginUser; 29 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 30 31 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 32 33 PetscCall(VecCreate(PETSC_COMM_WORLD, &f)); 34 PetscCall(VecSetSizes(f, 1, PETSC_DECIDE)); 35 PetscCall(VecSetFromOptions(f)); 36 PetscCall(VecSetUp(f)); 37 PetscCall(TSSetRHSFunction(ts, f, RHSFunction, NULL)); 38 PetscCall(VecDestroy(&f)); 39 40 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 41 PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE)); 42 PetscCall(MatSetFromOptions(A)); 43 PetscCall(MatSetUp(A)); 44 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 45 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 46 /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */ 47 PetscCall(MatShift(A, (PetscReal)1)); 48 PetscCall(MatShift(A, (PetscReal)-1)); 49 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL)); 50 PetscCall(MatDestroy(&A)); 51 52 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 53 PetscCall(VecSetSizes(x, 1, PETSC_DECIDE)); 54 PetscCall(VecSetFromOptions(x)); 55 PetscCall(VecSetUp(x)); 56 PetscCall(TSSetSolution(ts, x)); 57 PetscCall(VecDestroy(&x)); 58 59 PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL)); 60 PetscCall(TSSetPreStep(ts, PreStep)); 61 PetscCall(TSSetPostStep(ts, PostStep)); 62 63 { 64 TSAdapt adapt; 65 PetscCall(TSGetAdapt(ts, &adapt)); 66 PetscCall(TSAdaptSetType(adapt, TSADAPTNONE)); 67 } 68 { 69 PetscInt direction[3]; 70 PetscBool terminate[3]; 71 direction[0] = +1; 72 terminate[0] = PETSC_FALSE; 73 direction[1] = -1; 74 terminate[1] = PETSC_FALSE; 75 direction[2] = 0; 76 terminate[2] = PETSC_FALSE; 77 PetscCall(TSSetTimeStep(ts, 1)); 78 PetscCall(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL)); 79 } 80 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 81 PetscCall(TSSetResize(ts, PETSC_TRUE, TransferSetUp, Transfer, ntransfer)); 82 PetscCall(TSSetFromOptions(ts)); 83 84 /* --- First Solve --- */ 85 86 PetscCall(TSSetStepNumber(ts, 0)); 87 PetscCall(TSSetTimeStep(ts, 1)); 88 PetscCall(TSSetTime(ts, 0)); 89 PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL)); 90 PetscCall(TSSetMaxSteps(ts, 3)); 91 92 PetscCall(TSGetTime(ts, &t)); 93 PetscCall(TSGetSolution(ts, &x)); 94 PetscCall(VecSet(x, t)); 95 while (t < t_end) { 96 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n")); 97 PetscCall(TSSolve(ts, NULL)); 98 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n")); 99 PetscCall(TSGetTime(ts, &t)); 100 PetscCall(TSGetStepNumber(ts, &n)); 101 PetscCall(TSSetMaxSteps(ts, PetscMin(n + 3, n_end))); 102 } 103 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n")); 104 PetscCall(TSSolve(ts, NULL)); 105 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n")); 106 107 /* --- Second Solve --- */ 108 109 PetscCall(TSSetStepNumber(ts, 0)); 110 PetscCall(TSSetTimeStep(ts, 1)); 111 PetscCall(TSSetTime(ts, 0)); 112 PetscCall(TSSetMaxTime(ts, 3)); 113 PetscCall(TSSetMaxSteps(ts, PETSC_INT_MAX)); 114 115 PetscCall(TSGetTime(ts, &t)); 116 PetscCall(TSGetSolution(ts, &x)); 117 PetscCall(VecSet(x, t)); 118 while (t < t_end) { 119 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n")); 120 PetscCall(TSSolve(ts, NULL)); 121 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n")); 122 PetscCall(TSGetTime(ts, &t)); 123 PetscCall(TSSetMaxTime(ts, PetscMin(t + 3, t_end))); 124 } 125 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n")); 126 PetscCall(TSSolve(ts, NULL)); 127 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n")); 128 129 /* --- */ 130 131 PetscCall(TSDestroy(&ts)); 132 133 PetscCall(PetscFinalize()); 134 return 0; 135 } 136 137 /* -------------------------------------------------------------------*/ 138 139 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, PetscCtx ctx) 140 { 141 PetscFunctionBeginUser; 142 PetscCall(VecSet(f, (PetscReal)1)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, PetscCtx ctx) 147 { 148 PetscFunctionBeginUser; 149 PetscCall(MatZeroEntries(B)); 150 if (B != A) PetscCall(MatZeroEntries(A)); 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 PetscErrorCode PreStep(TS ts) 155 { 156 PetscInt n; 157 PetscReal t; 158 Vec x; 159 const PetscScalar *a; 160 PetscBool flg; 161 162 PetscFunctionBeginUser; 163 PetscCall(TSGetStepNumber(ts, &n)); 164 PetscCall(TSGetTime(ts, &t)); 165 PetscCall(TSGetSolution(ts, &x)); 166 PetscCall(VecGetArrayRead(x, &a)); 167 PetscCall(TSGetStepResize(ts, &flg)); 168 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g%s\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]), flg ? " resized" : "")); 169 PetscCall(VecRestoreArrayRead(x, &a)); 170 PetscFunctionReturn(PETSC_SUCCESS); 171 } 172 173 PetscErrorCode PostStep(TS ts) 174 { 175 PetscInt n; 176 PetscReal t; 177 Vec x; 178 const PetscScalar *a; 179 180 PetscFunctionBeginUser; 181 PetscCall(TSGetStepNumber(ts, &n)); 182 PetscCall(TSGetTime(ts, &t)); 183 PetscCall(TSGetSolution(ts, &x)); 184 PetscCall(VecGetArrayRead(x, &a)); 185 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]))); 186 PetscCall(VecRestoreArrayRead(x, &a)); 187 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189 190 PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, PetscCtx ctx) 191 { 192 const PetscScalar *a; 193 194 PetscFunctionBeginUser; 195 PetscCall(VecGetArrayRead(x, &a)); 196 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]))); 197 PetscCall(VecRestoreArrayRead(x, &a)); 198 PetscFunctionReturn(PETSC_SUCCESS); 199 } 200 201 PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscReal *fvalue, PetscCtx ctx) 202 { 203 PetscFunctionBeginUser; 204 fvalue[0] = t - 5; 205 fvalue[1] = 7 - t; 206 fvalue[2] = t - 9; 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, PetscCtx ctx) 211 { 212 PetscInt i; 213 const PetscScalar *a; 214 215 PetscFunctionBeginUser; 216 PetscCall(TSGetStepNumber(ts, &i)); 217 PetscCall(VecGetArrayRead(x, &a)); 218 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0]))); 219 PetscCall(VecRestoreArrayRead(x, &a)); 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 223 PetscErrorCode TransferSetUp(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, PetscCtx ctx) 224 { 225 PetscInt *nt = (PetscInt *)ctx; 226 227 PetscFunctionBeginUser; 228 if (n == 1) { 229 nt[0] = 2; 230 nt[1] = 2; 231 } 232 *flg = (PetscBool)(nt[0] && n && n % (nt[0]) == 0); 233 if (*flg) nt[0] += nt[1]; 234 PetscFunctionReturn(PETSC_SUCCESS); 235 } 236 237 PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], PetscCtx ctx) 238 { 239 PetscInt i; 240 241 PetscFunctionBeginUser; 242 PetscCall(TSGetStepNumber(ts, &i)); 243 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " nv %" PetscInt_FMT "\n", PETSC_FUNCTION_NAME, i, nv)); 244 for (i = 0; i < nv; i++) { 245 PetscCall(VecDuplicate(vin[i], &vout[i])); 246 PetscCall(VecCopy(vin[i], vout[i])); 247 } 248 PetscFunctionReturn(PETSC_SUCCESS); 249 } 250 251 /*TEST 252 253 test: 254 suffix: euler 255 diff_args: -j 256 args: -ts_type euler 257 output_file: output/ex1.out 258 259 test: 260 suffix: ssp 261 diff_args: -j 262 args: -ts_type ssp 263 output_file: output/ex1.out 264 265 test: 266 suffix: rk 267 diff_args: -j 268 args: -ts_type rk 269 output_file: output/ex1.out 270 271 test: 272 suffix: beuler 273 diff_args: -j 274 args: -ts_type beuler 275 output_file: output/ex1_theta.out 276 277 test: 278 suffix: cn 279 diff_args: -j 280 args: -ts_type cn 281 output_file: output/ex1_theta.out 282 283 test: 284 suffix: theta 285 args: -ts_type theta 286 diff_args: -j 287 output_file: output/ex1_theta.out 288 289 test: 290 suffix: bdf 291 diff_args: -j 292 args: -ts_type bdf 293 output_file: output/ex1_bdf.out 294 295 test: 296 suffix: alpha 297 diff_args: -j 298 args: -ts_type alpha 299 output_file: output/ex1_alpha.out 300 301 test: 302 suffix: rosw 303 diff_args: -j 304 args: -ts_type rosw 305 output_file: output/ex1.out 306 307 test: 308 suffix: arkimex 309 diff_args: -j 310 args: -ts_type arkimex 311 output_file: output/ex1_arkimex.out 312 313 TEST*/ 314