Actual source code: ex1.c


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

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

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

 10: static PetscErrorCode PreStep(TS);
 11: static PetscErrorCode PostStep(TS);
 12: static PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
 13: static PetscErrorCode Event(TS, PetscReal, Vec, PetscScalar *, void *);
 14: static PetscErrorCode PostEvent(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *);

 16: int main(int argc, char **argv)
 17: {
 18:   TS              ts;
 19:   PetscInt        n;
 20:   const PetscInt  n_end = 11;
 21:   PetscReal       t;
 22:   const PetscReal t_end = 11;
 23:   Vec             x;
 24:   Vec             f;
 25:   Mat             A;

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

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

 32:   PetscCall(VecCreate(PETSC_COMM_WORLD, &f));
 33:   PetscCall(VecSetSizes(f, 1, PETSC_DECIDE));
 34:   PetscCall(VecSetFromOptions(f));
 35:   PetscCall(VecSetUp(f));
 36:   PetscCall(TSSetRHSFunction(ts, f, RHSFunction, NULL));
 37:   PetscCall(VecDestroy(&f));

 39:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 40:   PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE));
 41:   PetscCall(MatSetFromOptions(A));
 42:   PetscCall(MatSetUp(A));
 43:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 44:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 45:   /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
 46:   PetscCall(MatShift(A, (PetscReal)1));
 47:   PetscCall(MatShift(A, (PetscReal)-1));
 48:   PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
 49:   PetscCall(MatDestroy(&A));

 51:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 52:   PetscCall(VecSetSizes(x, 1, PETSC_DECIDE));
 53:   PetscCall(VecSetFromOptions(x));
 54:   PetscCall(VecSetUp(x));
 55:   PetscCall(TSSetSolution(ts, x));
 56:   PetscCall(VecDestroy(&x));

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

 62:   {
 63:     TSAdapt adapt;
 64:     PetscCall(TSGetAdapt(ts, &adapt));
 65:     PetscCall(TSAdaptSetType(adapt, TSADAPTNONE));
 66:   }
 67:   {
 68:     PetscInt  direction[3];
 69:     PetscBool terminate[3];
 70:     direction[0] = +1;
 71:     terminate[0] = PETSC_FALSE;
 72:     direction[1] = -1;
 73:     terminate[1] = PETSC_FALSE;
 74:     direction[2] = 0;
 75:     terminate[2] = PETSC_FALSE;
 76:     PetscCall(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL));
 77:   }
 78:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
 79:   PetscCall(TSSetFromOptions(ts));

 81:   /* --- First Solve --- */

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

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

104:   /* --- Second Solve --- */

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

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

126:   /* --- */

128:   PetscCall(TSDestroy(&ts));

130:   PetscCall(PetscFinalize());
131:   return 0;
132: }

134: /* -------------------------------------------------------------------*/

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

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

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

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

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

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

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

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

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

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

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

218: /*TEST

220:     test:
221:       suffix: euler
222:       args: -ts_type euler
223:       output_file: output/ex1.out

225:     test:
226:       suffix: ssp
227:       args:   -ts_type ssp
228:       output_file: output/ex1.out

230:     test:
231:       suffix: rk
232:       args: -ts_type rk
233:       output_file: output/ex1.out

235:     test:
236:       suffix: beuler
237:       args: -ts_type beuler
238:       output_file: output/ex1.out

240:     test:
241:       suffix: cn
242:       args: -ts_type cn
243:       output_file: output/ex1.out

245:     test:
246:       suffix: theta
247:       args: -ts_type theta
248:       output_file: output/ex1.out

250:     test:
251:       suffix: bdf
252:       args: -ts_type bdf
253:       output_file: output/ex1.out

255:     test:
256:       suffix: alpha
257:       args: -ts_type alpha
258:       output_file: output/ex1.out

260:     test:
261:       suffix: rosw
262:       args: -ts_type rosw
263:       output_file: output/ex1.out

265:     test:
266:       suffix: arkimex
267:       args: -ts_type arkimex
268:       output_file: output/ex1.out

270: TEST*/