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; 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, 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 *flg = (PetscBool)(*nt && n && n % (*nt) == 0); 227 PetscFunctionReturn(PETSC_SUCCESS); 228 } 229 230 PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *ctx) 231 { 232 PetscInt i; 233 234 PetscFunctionBeginUser; 235 PetscCall(TSGetStepNumber(ts, &i)); 236 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " nv %" PetscInt_FMT "\n", PETSC_FUNCTION_NAME, i, nv)); 237 for (i = 0; i < nv; i++) { 238 PetscCall(VecDuplicate(vin[i], &vout[i])); 239 PetscCall(VecCopy(vin[i], vout[i])); 240 } 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 /*TEST 245 246 test: 247 suffix: euler 248 diff_args: -j 249 args: -ts_type euler 250 output_file: output/ex1.out 251 252 test: 253 suffix: ssp 254 diff_args: -j 255 args: -ts_type ssp 256 output_file: output/ex1.out 257 258 test: 259 suffix: rk 260 diff_args: -j 261 args: -ts_type rk 262 output_file: output/ex1.out 263 264 test: 265 suffix: beuler 266 diff_args: -j 267 args: -ts_type beuler 268 output_file: output/ex1.out 269 270 test: 271 suffix: cn 272 diff_args: -j 273 args: -ts_type cn 274 output_file: output/ex1.out 275 276 test: 277 suffix: theta 278 args: -ts_type theta 279 diff_args: -j 280 output_file: output/ex1.out 281 282 test: 283 suffix: bdf 284 diff_args: -j 285 args: -ts_type bdf 286 output_file: output/ex1_bdf.out 287 288 test: 289 suffix: alpha 290 diff_args: -j 291 args: -ts_type alpha 292 output_file: output/ex1.out 293 294 test: 295 suffix: rosw 296 diff_args: -j 297 args: -ts_type rosw 298 output_file: output/ex1.out 299 300 test: 301 suffix: arkimex 302 diff_args: -j 303 args: -ts_type arkimex 304 output_file: output/ex1.out 305 306 TEST*/ 307