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

 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;

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

 84:   /* --- First Solve --- */

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

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

107:   /* --- Second Solve --- */

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

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

129:   /* --- */

131:   PetscCall(TSDestroy(&ts));

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

137: /* -------------------------------------------------------------------*/

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: }

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: }

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

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: }

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

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: }

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

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: }

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: }

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;

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: }

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

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: }

235: PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *ctx)
236: {
237:   PetscInt i;

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: }

249: /*TEST

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

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

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

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

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

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

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

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

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

305:     test:
306:       suffix: arkimex
307:       diff_args: -j
308:       args: -ts_type arkimex
309:       output_file: output/ex1_arkimex.out

311: TEST*/