Actual source code: ex17.c
1: static char help[] = "Solves the ODE du/dt = poly(t), u(0) = 0. Tests TSResize for varying size.\n\n";
3: #include <petscts.h>
5: PetscScalar poly(PetscInt p, PetscReal t)
6: {
7: return p ? t * poly(p - 1, t) : 1.0;
8: }
10: PetscScalar dpoly(PetscInt p, PetscReal t)
11: {
12: return p > 0 ? (PetscReal)p * poly(p - 1, t) : 0.0;
13: }
15: PetscErrorCode CreateVec(PetscInt lsize, Vec *out)
16: {
17: Vec x;
19: PetscFunctionBeginUser;
20: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
21: PetscCall(VecSetSizes(x, lsize, PETSC_DECIDE));
22: PetscCall(VecSetFromOptions(x));
23: PetscCall(VecSetUp(x));
24: *out = x;
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: PetscErrorCode CreateMat(PetscInt lsize, Mat *out)
29: {
30: Mat A;
32: PetscFunctionBeginUser;
33: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
34: PetscCall(MatSetSizes(A, lsize, lsize, PETSC_DECIDE, PETSC_DECIDE));
35: PetscCall(MatSetFromOptions(A));
36: PetscCall(MatSetUp(A));
37: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
38: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
39: /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
40: PetscCall(MatShift(A, (PetscReal)1));
41: PetscCall(MatShift(A, (PetscReal)-1));
42: *out = A;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
47: {
48: PetscInt *order = (PetscInt *)ctx;
50: PetscFunctionBeginUser;
51: PetscCall(VecSet(f, dpoly(*order, t)));
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
56: {
57: PetscFunctionBeginUser;
58: PetscCall(MatZeroEntries(B));
59: if (B != A) PetscCall(MatZeroEntries(A));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
64: {
65: PetscInt n, nnew;
67: PetscFunctionBeginUser;
68: PetscAssert(nv > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Zero vectors");
69: PetscCall(VecGetLocalSize(vecsin[0], &n));
70: nnew = n == 2 ? 1 : 2;
71: for (PetscInt i = 0; i < nv; i++) {
72: const PetscScalar *vals;
74: PetscCall(CreateVec(nnew, &vecsout[i]));
75: PetscCall(VecGetArrayRead(vecsin[i], &vals));
76: PetscCall(VecSet(vecsout[i], vals[0]));
77: PetscCall(VecRestoreArrayRead(vecsin[i], &vals));
78: }
79: Mat A;
80: PetscCall(CreateMat(nnew, &A));
81: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
82: PetscCall(MatDestroy(&A));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PetscErrorCode TransferSetUp(TS ts, PetscInt step, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
87: {
88: PetscFunctionBeginUser;
89: *resize = PETSC_TRUE;
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
94: {
95: const PetscScalar *a;
96: PetscScalar *store = (PetscScalar *)ctx;
98: PetscFunctionBeginUser;
99: PetscCall(VecGetArrayRead(x, &a));
100: if (n < 10) store[n] = a[0];
101: PetscCall(VecRestoreArrayRead(x, &a));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: int main(int argc, char **argv)
106: {
107: TS ts;
108: Vec x;
109: Mat A;
110: PetscInt order = 2;
111: PetscScalar results[2][10];
112: /* I would like to use 0 here, but linux-gcc-complex-opt-32bit errors with arkimex with 1.e-18 errors, macOS clang requires an even larger tolerance */
113: PetscReal tol = 10 * PETSC_MACHINE_EPSILON;
115: PetscFunctionBeginUser;
116: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
117: PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL));
118: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));
120: for (PetscInt i = 0; i < 2; i++) {
121: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
122: PetscCall(TSSetProblemType(ts, TS_LINEAR));
124: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &order));
126: PetscCall(CreateMat(1, &A));
127: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
128: PetscCall(MatDestroy(&A));
130: PetscCall(CreateVec(1, &x));
131: PetscCall(TSSetSolution(ts, x));
132: PetscCall(VecDestroy(&x));
134: for (PetscInt j = 0; j < 10; j++) results[i][j] = 0;
135: PetscCall(TSMonitorSet(ts, Monitor, results[i], NULL));
136: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
137: if (i) PetscCall(TSSetResize(ts, TransferSetUp, Transfer, NULL));
138: PetscCall(TSSetTime(ts, 0));
139: PetscCall(TSSetTimeStep(ts, 1. / 4.));
140: PetscCall(TSSetMaxSteps(ts, 10));
141: PetscCall(TSSetFromOptions(ts));
143: PetscCall(TSSolve(ts, NULL));
145: PetscCall(TSDestroy(&ts));
146: }
148: /* Dump errors if any */
149: PetscBool flg = PETSC_FALSE;
150: for (PetscInt i = 0; i < 10; i++) {
151: PetscReal err = PetscAbsScalar(results[0][i] - results[1][i]);
152: if (err > tol) {
153: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error step %" PetscInt_FMT ": %g\n", i, (double)err));
154: flg = PETSC_TRUE;
155: }
156: }
157: if (flg) {
158: PetscCall(PetscScalarView(10, results[0], PETSC_VIEWER_STDOUT_WORLD));
159: PetscCall(PetscScalarView(10, results[1], PETSC_VIEWER_STDOUT_WORLD));
160: }
161: PetscCall(PetscFinalize());
162: return 0;
163: }
165: /*TEST
167: test:
168: suffix: bdf
169: args: -ts_adapt_wnormtype infinity -ts_type bdf -ts_bdf_order {{2 3 4 5 6}} -order 6 -ts_adapt_type {{none basic dsp}} -ksp_type preonly -pc_type lu
170: output_file: output/ex17.out
172: test:
173: suffix: expl
174: args: -ts_adapt_wnormtype infinity -ts_type {{euler rk ssp}} -order 6 -ts_adapt_type {{none basic dsp}}
175: output_file: output/ex17.out
177: test:
178: suffix: impl
179: args: -ts_adapt_wnormtype infinity -ts_type {{rosw beuler cn alpha theta arkimex}} -order 6 -ts_adapt_type {{none basic dsp}} -ksp_type preonly -pc_type lu
180: output_file: output/ex17.out
182: TEST*/