1 static char help[] = "Solves the ODE du/dt = poly(t), u(0) = 0. Tests TSResize for varying size.\n\n";
2
3 #include <petscts.h>
4
poly(PetscInt p,PetscReal t)5 PetscScalar poly(PetscInt p, PetscReal t)
6 {
7 return p ? t * poly(p - 1, t) : 1.0;
8 }
9
dpoly(PetscInt p,PetscReal t)10 PetscScalar dpoly(PetscInt p, PetscReal t)
11 {
12 return p > 0 ? (PetscReal)p * poly(p - 1, t) : 0.0;
13 }
14
CreateVec(PetscInt lsize,Vec * out)15 PetscErrorCode CreateVec(PetscInt lsize, Vec *out)
16 {
17 Vec x;
18
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 }
27
CreateMat(PetscInt lsize,Mat * out)28 PetscErrorCode CreateMat(PetscInt lsize, Mat *out)
29 {
30 Mat A;
31
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 }
45
RHSFunction(TS ts,PetscReal t,Vec x,Vec f,PetscCtx ctx)46 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, PetscCtx ctx)
47 {
48 PetscInt *order = (PetscInt *)ctx;
49
50 PetscFunctionBeginUser;
51 PetscCall(VecSet(f, dpoly(*order, t)));
52 PetscFunctionReturn(PETSC_SUCCESS);
53 }
54
RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat B,PetscCtx ctx)55 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, PetscCtx ctx)
56 {
57 PetscFunctionBeginUser;
58 PetscCall(MatZeroEntries(B));
59 if (B != A) PetscCall(MatZeroEntries(A));
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
Transfer(TS ts,PetscInt nv,Vec vecsin[],Vec vecsout[],PetscCtx ctx)63 PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], PetscCtx ctx)
64 {
65 PetscInt n, nnew;
66
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;
73
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 }
85
TransferSetUp(TS ts,PetscInt step,PetscReal time,Vec sol,PetscBool * resize,PetscCtx ctx)86 PetscErrorCode TransferSetUp(TS ts, PetscInt step, PetscReal time, Vec sol, PetscBool *resize, PetscCtx ctx)
87 {
88 PetscBool *alreadydone = (PetscBool *)ctx;
89
90 PetscFunctionBeginUser;
91 *alreadydone = (PetscBool)!*alreadydone;
92 *resize = *alreadydone;
93 PetscFunctionReturn(PETSC_SUCCESS);
94 }
95
Monitor(TS ts,PetscInt n,PetscReal t,Vec x,PetscCtx ctx)96 PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, PetscCtx ctx)
97 {
98 const PetscScalar *a;
99 PetscScalar *store = (PetscScalar *)ctx;
100
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 }
107
main(int argc,char ** argv)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;
117
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));
122
123 for (PetscInt i = 0; i < 3; i++) {
124 PetscBool alreadydone = PETSC_TRUE;
125
126 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
127 PetscCall(TSSetProblemType(ts, TS_LINEAR));
128
129 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &order));
130
131 PetscCall(CreateMat(1, &A));
132 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
133 PetscCall(MatDestroy(&A));
134
135 PetscCall(CreateVec(1, &x));
136 PetscCall(TSSetSolution(ts, x));
137 PetscCall(VecDestroy(&x));
138
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));
147
148 PetscCall(TSSolve(ts, NULL));
149
150 PetscCall(TSDestroy(&ts));
151 }
152
153 /* Dump errors if any */
154 PetscBool flg = PETSC_FALSE;
155 for (PetscInt i = 0; i < 10; i++) {
156 PetscReal err;
157
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 }
177
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
182
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/empty.out
187
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/empty.out
192
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/empty.out
197
198 TEST*/
199