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, (char *)0, 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_MAX_INT)); 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, void *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, void *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 161 PetscFunctionBeginUser; 162 PetscCall(TSGetStepNumber(ts, &n)); 163 PetscCall(TSGetTime(ts, &t)); 164 PetscCall(TSGetSolution(ts, &x)); 165 PetscCall(VecGetArrayRead(x, &a)); 166 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]))); 167 PetscCall(VecRestoreArrayRead(x, &a)); 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170 171 PetscErrorCode PostStep(TS ts) 172 { 173 PetscInt n; 174 PetscReal t; 175 Vec x; 176 const PetscScalar *a; 177 178 PetscFunctionBeginUser; 179 PetscCall(TSGetStepNumber(ts, &n)); 180 PetscCall(TSGetTime(ts, &t)); 181 PetscCall(TSGetSolution(ts, &x)); 182 PetscCall(VecGetArrayRead(x, &a)); 183 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]))); 184 PetscCall(VecRestoreArrayRead(x, &a)); 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187 188 PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx) 189 { 190 const PetscScalar *a; 191 192 PetscFunctionBeginUser; 193 PetscCall(VecGetArrayRead(x, &a)); 194 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]))); 195 PetscCall(VecRestoreArrayRead(x, &a)); 196 PetscFunctionReturn(PETSC_SUCCESS); 197 } 198 199 PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscReal *fvalue, void *ctx) 200 { 201 PetscFunctionBeginUser; 202 fvalue[0] = t - 5; 203 fvalue[1] = 7 - t; 204 fvalue[2] = t - 9; 205 PetscFunctionReturn(PETSC_SUCCESS); 206 } 207 208 PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx) 209 { 210 PetscInt i; 211 const PetscScalar *a; 212 213 PetscFunctionBeginUser; 214 PetscCall(TSGetStepNumber(ts, &i)); 215 PetscCall(VecGetArrayRead(x, &a)); 216 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0]))); 217 PetscCall(VecRestoreArrayRead(x, &a)); 218 PetscFunctionReturn(PETSC_SUCCESS); 219 } 220 221 PetscErrorCode TransferSetUp(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *ctx) 222 { 223 PetscInt *nt = (PetscInt *)ctx; 224 225 PetscFunctionBeginUser; 226 if (n == 1) { 227 nt[0] = 2; 228 nt[1] = 2; 229 } 230 *flg = (PetscBool)(nt[0] && n && n % (nt[0]) == 0); 231 if (*flg) nt[0] += nt[1]; 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *ctx) 236 { 237 PetscInt i; 238 239 PetscFunctionBeginUser; 240 PetscCall(TSGetStepNumber(ts, &i)); 241 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " nv %" PetscInt_FMT "\n", PETSC_FUNCTION_NAME, i, nv)); 242 for (i = 0; i < nv; i++) { 243 PetscCall(VecDuplicate(vin[i], &vout[i])); 244 PetscCall(VecCopy(vin[i], vout[i])); 245 } 246 PetscFunctionReturn(PETSC_SUCCESS); 247 } 248 249 /*TEST 250 251 test: 252 suffix: euler 253 diff_args: -j 254 args: -ts_type euler 255 output_file: output/ex1.out 256 257 test: 258 suffix: ssp 259 diff_args: -j 260 args: -ts_type ssp 261 output_file: output/ex1.out 262 263 test: 264 suffix: rk 265 diff_args: -j 266 args: -ts_type rk 267 output_file: output/ex1.out 268 269 test: 270 suffix: beuler 271 diff_args: -j 272 args: -ts_type beuler 273 output_file: output/ex1_theta.out 274 275 test: 276 suffix: cn 277 diff_args: -j 278 args: -ts_type cn 279 output_file: output/ex1_theta.out 280 281 test: 282 suffix: theta 283 args: -ts_type theta 284 diff_args: -j 285 output_file: output/ex1_theta.out 286 287 test: 288 suffix: bdf 289 diff_args: -j 290 args: -ts_type bdf 291 output_file: output/ex1_bdf.out 292 293 test: 294 suffix: alpha 295 diff_args: -j 296 args: -ts_type alpha 297 output_file: output/ex1_alpha.out 298 299 test: 300 suffix: rosw 301 diff_args: -j 302 args: -ts_type rosw 303 output_file: output/ex1.out 304 305 test: 306 suffix: arkimex 307 diff_args: -j 308 args: -ts_type arkimex 309 output_file: output/ex1_arkimex.out 310 311 TEST*/ 312