Actual source code: ex1.c

  1: static char help[] = "Solves the trivial ODE du/dt = 1, u(0) = 0. \n\n";

  3: #include <petscts.h>
  4: #include <petscpc.h>

  6: static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
  7: static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);

  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, PetscScalar *, 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 *);

 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;

 28:   PetscFunctionBeginUser;
 29:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 31:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));

 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));

 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));

 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));

 59:   PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
 60:   PetscCall(TSSetPreStep(ts, PreStep));
 61:   PetscCall(TSSetPostStep(ts, PostStep));

 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(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL));
 78:   }
 79:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
 80:   PetscCall(TSSetResize(ts, TransferSetUp, Transfer, &ntransfer));
 81:   PetscCall(TSSetFromOptions(ts));

 83:   /* --- First Solve --- */

 85:   PetscCall(TSSetStepNumber(ts, 0));
 86:   PetscCall(TSSetTimeStep(ts, 1));
 87:   PetscCall(TSSetTime(ts, 0));
 88:   PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL));
 89:   PetscCall(TSSetMaxSteps(ts, 3));

 91:   PetscCall(TSGetTime(ts, &t));
 92:   PetscCall(TSGetSolution(ts, &x));
 93:   PetscCall(VecSet(x, t));
 94:   while (t < t_end) {
 95:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
 96:     PetscCall(TSSolve(ts, NULL));
 97:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
 98:     PetscCall(TSGetTime(ts, &t));
 99:     PetscCall(TSGetStepNumber(ts, &n));
100:     PetscCall(TSSetMaxSteps(ts, PetscMin(n + 3, n_end)));
101:   }
102:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
103:   PetscCall(TSSolve(ts, NULL));
104:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));

106:   /* --- Second Solve --- */

108:   PetscCall(TSSetStepNumber(ts, 0));
109:   PetscCall(TSSetTimeStep(ts, 1));
110:   PetscCall(TSSetTime(ts, 0));
111:   PetscCall(TSSetMaxTime(ts, 3));
112:   PetscCall(TSSetMaxSteps(ts, PETSC_MAX_INT));

114:   PetscCall(TSGetTime(ts, &t));
115:   PetscCall(TSGetSolution(ts, &x));
116:   PetscCall(VecSet(x, t));
117:   while (t < t_end) {
118:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
119:     PetscCall(TSSolve(ts, NULL));
120:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
121:     PetscCall(TSGetTime(ts, &t));
122:     PetscCall(TSSetMaxTime(ts, PetscMin(t + 3, t_end)));
123:   }
124:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
125:   PetscCall(TSSolve(ts, NULL));
126:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));

128:   /* --- */

130:   PetscCall(TSDestroy(&ts));

132:   PetscCall(PetscFinalize());
133:   return 0;
134: }

136: /* -------------------------------------------------------------------*/

138: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
139: {
140:   PetscFunctionBeginUser;
141:   PetscCall(VecSet(f, (PetscReal)1));
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }

145: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
146: {
147:   PetscFunctionBeginUser;
148:   PetscCall(MatZeroEntries(B));
149:   if (B != A) PetscCall(MatZeroEntries(A));
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: PetscErrorCode PreStep(TS ts)
154: {
155:   PetscInt           n;
156:   PetscReal          t;
157:   Vec                x;
158:   const PetscScalar *a;

160:   PetscFunctionBeginUser;
161:   PetscCall(TSGetStepNumber(ts, &n));
162:   PetscCall(TSGetTime(ts, &t));
163:   PetscCall(TSGetSolution(ts, &x));
164:   PetscCall(VecGetArrayRead(x, &a));
165:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
166:   PetscCall(VecRestoreArrayRead(x, &a));
167:   PetscFunctionReturn(PETSC_SUCCESS);
168: }

170: PetscErrorCode PostStep(TS ts)
171: {
172:   PetscInt           n;
173:   PetscReal          t;
174:   Vec                x;
175:   const PetscScalar *a;

177:   PetscFunctionBeginUser;
178:   PetscCall(TSGetStepNumber(ts, &n));
179:   PetscCall(TSGetTime(ts, &t));
180:   PetscCall(TSGetSolution(ts, &x));
181:   PetscCall(VecGetArrayRead(x, &a));
182:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
183:   PetscCall(VecRestoreArrayRead(x, &a));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
188: {
189:   const PetscScalar *a;

191:   PetscFunctionBeginUser;
192:   PetscCall(VecGetArrayRead(x, &a));
193:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
194:   PetscCall(VecRestoreArrayRead(x, &a));
195:   PetscFunctionReturn(PETSC_SUCCESS);
196: }

198: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
199: {
200:   PetscFunctionBeginUser;
201:   fvalue[0] = t - 5;
202:   fvalue[1] = 7 - t;
203:   fvalue[2] = t - 9;
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }

207: PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx)
208: {
209:   PetscInt           i;
210:   const PetscScalar *a;

212:   PetscFunctionBeginUser;
213:   PetscCall(TSGetStepNumber(ts, &i));
214:   PetscCall(VecGetArrayRead(x, &a));
215:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0])));
216:   PetscCall(VecRestoreArrayRead(x, &a));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: PetscErrorCode TransferSetUp(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *ctx)
221: {
222:   PetscInt *nt = (PetscInt *)ctx;

224:   PetscFunctionBeginUser;
225:   *flg = (PetscBool)(*nt && n && n % (*nt) == 0);
226:   PetscFunctionReturn(PETSC_SUCCESS);
227: }

229: PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *ctx)
230: {
231:   PetscInt i;

233:   PetscFunctionBeginUser;
234:   PetscCall(TSGetStepNumber(ts, &i));
235:   PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " nv %" PetscInt_FMT "\n", PETSC_FUNCTION_NAME, i, nv));
236:   for (i = 0; i < nv; i++) {
237:     PetscCall(VecDuplicate(vin[i], &vout[i]));
238:     PetscCall(VecCopy(vin[i], vout[i]));
239:   }
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: /*TEST

245:     test:
246:       suffix: euler
247:       diff_args: -j
248:       args: -ts_type euler
249:       output_file: output/ex1.out

251:     test:
252:       suffix: ssp
253:       diff_args: -j
254:       args: -ts_type ssp
255:       output_file: output/ex1.out

257:     test:
258:       suffix: rk
259:       diff_args: -j
260:       args: -ts_type rk
261:       output_file: output/ex1.out

263:     test:
264:       suffix: beuler
265:       diff_args: -j
266:       args: -ts_type beuler
267:       output_file: output/ex1.out

269:     test:
270:       suffix: cn
271:       diff_args: -j
272:       args: -ts_type cn
273:       output_file: output/ex1.out

275:     test:
276:       suffix: theta
277:       args: -ts_type theta
278:       diff_args: -j
279:       output_file: output/ex1.out

281:     test:
282:       suffix: bdf
283:       diff_args: -j
284:       args: -ts_type bdf
285:       output_file: output/ex1_bdf.out

287:     test:
288:       suffix: alpha
289:       diff_args: -j
290:       args: -ts_type alpha
291:       output_file: output/ex1.out

293:     test:
294:       suffix: rosw
295:       diff_args: -j
296:       args: -ts_type rosw
297:       output_file: output/ex1.out

299:     test:
300:       suffix: arkimex
301:       diff_args: -j
302:       args: -ts_type arkimex
303:       output_file: output/ex1.out

305: TEST*/