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;

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

 30:   TSCreate(PETSC_COMM_WORLD, &ts);

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

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

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

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

 62:   {
 63:     TSAdapt adapt;
 64:     TSGetAdapt(ts, &adapt);
 65:     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:     TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL);
 77:   }
 78:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
 79:   TSSetFromOptions(ts);

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

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

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

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

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

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

126:   /* --- */

128:   TSDestroy(&ts);

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

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

136: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
137: {
139:   VecSet(f, (PetscReal)1);
140:   return 0;
141: }

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

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

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

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

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

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

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

196: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
197: {
199:   fvalue[0] = t - 5;
200:   fvalue[1] = 7 - t;
201:   fvalue[2] = t - 9;
202:   return 0;
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;

211:   TSGetStepNumber(ts, &i);
212:   VecGetArrayRead(x, &a);
213:   PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0]));
214:   VecRestoreArrayRead(x, &a);
215:   return 0;
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*/