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: PetscBool *alreadydone = (PetscBool *)ctx;
90: PetscFunctionBeginUser;
91: *alreadydone = (PetscBool)!*alreadydone;
92: *resize = *alreadydone;
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
97: {
98: const PetscScalar *a;
99: PetscScalar *store = (PetscScalar *)ctx;
101: PetscFunctionBeginUser;
102: PetscCall(VecGetArrayRead(x, &a));
103: if (n < 10) store[n] = a[0];
104: PetscCall(VecRestoreArrayRead(x, &a));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: int main(int argc, char **argv)
109: {
110: TS ts;
111: Vec x;
112: Mat A;
113: PetscInt order = 2;
114: PetscScalar results[3][10];
115: /* I would like to use 0 here, but arkimex errors with 1.e-14 discrepancy when using TSRResize without restart on some machines (mostly arm-osx) */
116: PetscReal tol = PETSC_SMALL;
118: PetscFunctionBeginUser;
119: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
120: PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL));
121: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));
123: for (PetscInt i = 0; i < 3; i++) {
124: PetscBool alreadydone = PETSC_TRUE;
126: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
127: PetscCall(TSSetProblemType(ts, TS_LINEAR));
129: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &order));
131: PetscCall(CreateMat(1, &A));
132: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
133: PetscCall(MatDestroy(&A));
135: PetscCall(CreateVec(1, &x));
136: PetscCall(TSSetSolution(ts, x));
137: PetscCall(VecDestroy(&x));
139: for (PetscInt j = 0; j < 10; j++) results[i][j] = 0;
140: PetscCall(TSMonitorSet(ts, Monitor, results[i], NULL));
141: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
142: if (i) PetscCall(TSSetResize(ts, i == 1 ? PETSC_TRUE : PETSC_FALSE, TransferSetUp, Transfer, &alreadydone));
143: PetscCall(TSSetTime(ts, 0));
144: PetscCall(TSSetTimeStep(ts, 1. / 4.));
145: PetscCall(TSSetMaxSteps(ts, 10));
146: PetscCall(TSSetFromOptions(ts));
148: PetscCall(TSSolve(ts, NULL));
150: PetscCall(TSDestroy(&ts));
151: }
153: /* Dump errors if any */
154: PetscBool flg = PETSC_FALSE;
155: for (PetscInt i = 0; i < 10; i++) {
156: PetscReal err;
158: err = PetscAbsScalar(results[0][i] - results[1][i]);
159: if (err > tol) {
160: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error with restart for step %" PetscInt_FMT ": %g\n", i, (double)err));
161: flg = PETSC_TRUE;
162: }
163: err = PetscAbsScalar(results[0][i] - results[2][i]);
164: if (err > tol) {
165: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error without restart for step %" PetscInt_FMT ": %g\n", i, (double)err));
166: flg = PETSC_TRUE;
167: }
168: }
169: if (flg) {
170: PetscCall(PetscScalarView(10, results[0], PETSC_VIEWER_STDOUT_WORLD));
171: PetscCall(PetscScalarView(10, results[1], PETSC_VIEWER_STDOUT_WORLD));
172: PetscCall(PetscScalarView(10, results[2], PETSC_VIEWER_STDOUT_WORLD));
173: }
174: PetscCall(PetscFinalize());
175: return 0;
176: }
178: /*TEST
179: testset:
180: # using gemv gives larger error which failed error checking
181: args: -vec_mdot_use_gemv 0 -vec_maxpy_use_gemv 0
183: test:
184: suffix: bdf
185: 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
186: output_file: output/ex17.out
188: test:
189: suffix: expl
190: args: -ts_adapt_wnormtype infinity -ts_type {{euler rk ssp}} -order 6 -ts_adapt_type {{none basic dsp}}
191: output_file: output/ex17.out
193: test:
194: suffix: impl
195: 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
196: output_file: output/ex17.out
198: TEST*/