Actual source code: ex2.c
1: /*
2: Formatted test for TS routines.
4: Solves U_t=F(t,u)
5: Where:
7: [2*u1+u2 ]
8: F(t,u)= [u1+2*u2+u3]
9: [ u2+2*u3]
10: We can compare the solutions from euler, beuler and SUNDIALS to
11: see what is the difference.
13: */
15: static char help[] = "Solves a linear ODE. \n\n";
17: #include <petscts.h>
18: #include <petscpc.h>
20: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
21: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
22: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
23: extern PetscErrorCode Initial(Vec, void *);
24: extern PetscErrorCode MyMatMult(Mat, Vec, Vec);
26: extern PetscReal solx(PetscReal);
27: extern PetscReal soly(PetscReal);
28: extern PetscReal solz(PetscReal);
30: int main(int argc, char **argv)
31: {
32: PetscInt time_steps = 100, steps;
33: Vec global;
34: PetscReal dt, ftime;
35: TS ts;
36: Mat A = 0, S;
37: PetscBool nest;
39: PetscFunctionBeginUser;
40: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
41: PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
42: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nest", &nest, NULL));
44: /* create vector to hold state */
45: if (nest) {
46: Vec g[3];
48: PetscCall(VecCreate(PETSC_COMM_WORLD, &g[0]));
49: PetscCall(VecSetSizes(g[0], PETSC_DECIDE, 1));
50: PetscCall(VecSetFromOptions(g[0]));
51: PetscCall(VecDuplicate(g[0], &g[1]));
52: PetscCall(VecDuplicate(g[0], &g[2]));
53: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 3, NULL, g, &global));
54: PetscCall(VecDestroy(&g[0]));
55: PetscCall(VecDestroy(&g[1]));
56: PetscCall(VecDestroy(&g[2]));
57: } else {
58: PetscCall(VecCreate(PETSC_COMM_WORLD, &global));
59: PetscCall(VecSetSizes(global, PETSC_DECIDE, 3));
60: PetscCall(VecSetFromOptions(global));
61: }
63: /* set initial conditions */
64: PetscCall(Initial(global, NULL));
66: /* make timestep context */
67: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
68: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
69: PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
70: dt = 0.001;
72: /*
73: The user provides the RHS and Jacobian
74: */
75: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
76: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
77: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3));
78: PetscCall(MatSetFromOptions(A));
79: PetscCall(MatSetUp(A));
80: PetscCall(RHSJacobian(ts, 0.0, global, A, A, NULL));
81: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
83: PetscCall(MatCreateShell(PETSC_COMM_WORLD, 3, 3, 3, 3, NULL, &S));
84: PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MyMatMult));
85: PetscCall(TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL));
87: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
88: PetscCall(TSSetFromOptions(ts));
90: PetscCall(TSSetTimeStep(ts, dt));
91: PetscCall(TSSetMaxSteps(ts, time_steps));
92: PetscCall(TSSetMaxTime(ts, 1));
93: PetscCall(TSSetSolution(ts, global));
95: PetscCall(TSSolve(ts, global));
96: PetscCall(TSGetSolveTime(ts, &ftime));
97: PetscCall(TSGetStepNumber(ts, &steps));
99: /* free the memories */
101: PetscCall(TSDestroy(&ts));
102: PetscCall(VecDestroy(&global));
103: PetscCall(MatDestroy(&A));
104: PetscCall(MatDestroy(&S));
106: PetscCall(PetscFinalize());
107: return 0;
108: }
110: PetscErrorCode MyMatMult(Mat S, Vec x, Vec y)
111: {
112: const PetscScalar *inptr;
113: PetscScalar *outptr;
115: PetscFunctionBeginUser;
116: PetscCall(VecGetArrayRead(x, &inptr));
117: PetscCall(VecGetArrayWrite(y, &outptr));
119: outptr[0] = 2.0 * inptr[0] + inptr[1];
120: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
121: outptr[2] = inptr[1] + 2.0 * inptr[2];
123: PetscCall(VecRestoreArrayRead(x, &inptr));
124: PetscCall(VecRestoreArrayWrite(y, &outptr));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /* this test problem has initial values (1,1,1). */
129: PetscErrorCode Initial(Vec global, void *ctx)
130: {
131: PetscScalar *localptr;
132: PetscInt i, mybase, myend, locsize;
134: PetscFunctionBeginUser;
135: /* determine starting point of each processor */
136: PetscCall(VecGetOwnershipRange(global, &mybase, &myend));
137: PetscCall(VecGetLocalSize(global, &locsize));
139: /* Initialize the array */
140: PetscCall(VecGetArrayWrite(global, &localptr));
141: for (i = 0; i < locsize; i++) localptr[i] = 1.0;
143: if (mybase == 0) localptr[0] = 1.0;
145: PetscCall(VecRestoreArrayWrite(global, &localptr));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx)
150: {
151: VecScatter scatter;
152: IS from, to;
153: PetscInt i, n, *idx;
154: Vec tmp_vec;
155: const PetscScalar *tmp;
157: PetscFunctionBeginUser;
158: /* Get the size of the vector */
159: PetscCall(VecGetSize(global, &n));
161: /* Set the index sets */
162: PetscCall(PetscMalloc1(n, &idx));
163: for (i = 0; i < n; i++) idx[i] = i;
165: /* Create local sequential vectors */
166: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec));
168: /* Create scatter context */
169: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from));
170: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to));
171: PetscCall(VecScatterCreate(global, from, tmp_vec, to, &scatter));
172: PetscCall(VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD));
173: PetscCall(VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD));
175: PetscCall(VecGetArrayRead(tmp_vec, &tmp));
176: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e u = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0]), (double)PetscRealPart(tmp[1]), (double)PetscRealPart(tmp[2])));
177: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e errors = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0] - solx(time)), (double)PetscRealPart(tmp[1] - soly(time)), (double)PetscRealPart(tmp[2] - solz(time))));
178: PetscCall(VecRestoreArrayRead(tmp_vec, &tmp));
179: PetscCall(VecScatterDestroy(&scatter));
180: PetscCall(ISDestroy(&from));
181: PetscCall(ISDestroy(&to));
182: PetscCall(PetscFree(idx));
183: PetscCall(VecDestroy(&tmp_vec));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
188: {
189: PetscScalar *outptr;
190: const PetscScalar *inptr;
191: PetscInt i, n, *idx;
192: IS from, to;
193: VecScatter scatter;
194: Vec tmp_in, tmp_out;
196: PetscFunctionBeginUser;
197: /* Get the length of parallel vector */
198: PetscCall(VecGetSize(globalin, &n));
200: /* Set the index sets */
201: PetscCall(PetscMalloc1(n, &idx));
202: for (i = 0; i < n; i++) idx[i] = i;
204: /* Create local sequential vectors */
205: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_in));
206: PetscCall(VecDuplicate(tmp_in, &tmp_out));
208: /* Create scatter context */
209: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from));
210: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to));
211: PetscCall(VecScatterCreate(globalin, from, tmp_in, to, &scatter));
212: PetscCall(VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD));
213: PetscCall(VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD));
214: PetscCall(VecScatterDestroy(&scatter));
216: /*Extract income array */
217: PetscCall(VecGetArrayRead(tmp_in, &inptr));
219: /* Extract outcome array*/
220: PetscCall(VecGetArrayWrite(tmp_out, &outptr));
222: outptr[0] = 2.0 * inptr[0] + inptr[1];
223: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
224: outptr[2] = inptr[1] + 2.0 * inptr[2];
226: PetscCall(VecRestoreArrayRead(tmp_in, &inptr));
227: PetscCall(VecRestoreArrayWrite(tmp_out, &outptr));
229: PetscCall(VecScatterCreate(tmp_out, from, globalout, to, &scatter));
230: PetscCall(VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));
231: PetscCall(VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));
233: /* Destroy idx and scatter */
234: PetscCall(ISDestroy(&from));
235: PetscCall(ISDestroy(&to));
236: PetscCall(VecScatterDestroy(&scatter));
237: PetscCall(VecDestroy(&tmp_in));
238: PetscCall(VecDestroy(&tmp_out));
239: PetscCall(PetscFree(idx));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx)
244: {
245: PetscScalar v[3];
246: const PetscScalar *tmp;
247: PetscInt idx[3], i;
249: PetscFunctionBeginUser;
250: idx[0] = 0;
251: idx[1] = 1;
252: idx[2] = 2;
253: PetscCall(VecGetArrayRead(x, &tmp));
255: i = 0;
256: v[0] = 2.0;
257: v[1] = 1.0;
258: v[2] = 0.0;
259: PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES));
261: i = 1;
262: v[0] = 1.0;
263: v[1] = 2.0;
264: v[2] = 1.0;
265: PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES));
267: i = 2;
268: v[0] = 0.0;
269: v[1] = 1.0;
270: v[2] = 2.0;
271: PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES));
273: PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
274: PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));
276: if (A != BB) {
277: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
278: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
279: }
280: PetscCall(VecRestoreArrayRead(x, &tmp));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*
285: The exact solutions
286: */
287: PetscReal solx(PetscReal t)
288: {
289: return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0));
290: }
292: PetscReal soly(PetscReal t)
293: {
294: return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0);
295: }
297: PetscReal solz(PetscReal t)
298: {
299: return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0));
300: }
302: /*TEST
304: test:
305: suffix: euler
306: args: -ts_type euler -nest {{0 1}}
307: requires: !single
309: test:
310: suffix: beuler
311: args: -ts_type beuler -nest {{0 1}}
312: requires: !single
314: test:
315: suffix: rk
316: args: -ts_type rk -nest {{0 1}} -ts_adapt_monitor
317: requires: !single
319: TEST*/