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;
39: PetscInitialize(&argc, &argv, (char *)0, help);
40: PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL);
42: /* set initial conditions */
43: VecCreate(PETSC_COMM_WORLD, &global);
44: VecSetSizes(global, PETSC_DECIDE, 3);
45: VecSetFromOptions(global);
46: Initial(global, NULL);
48: /* make timestep context */
49: TSCreate(PETSC_COMM_WORLD, &ts);
50: TSSetProblemType(ts, TS_NONLINEAR);
51: TSMonitorSet(ts, Monitor, NULL, NULL);
52: dt = 0.001;
54: /*
55: The user provides the RHS and Jacobian
56: */
57: TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
58: MatCreate(PETSC_COMM_WORLD, &A);
59: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3);
60: MatSetFromOptions(A);
61: MatSetUp(A);
62: RHSJacobian(ts, 0.0, global, A, A, NULL);
63: TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL);
65: MatCreateShell(PETSC_COMM_WORLD, 3, 3, 3, 3, NULL, &S);
66: MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MyMatMult);
67: TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL);
69: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
70: TSSetFromOptions(ts);
72: TSSetTimeStep(ts, dt);
73: TSSetMaxSteps(ts, time_steps);
74: TSSetMaxTime(ts, 1);
75: TSSetSolution(ts, global);
77: TSSolve(ts, global);
78: TSGetSolveTime(ts, &ftime);
79: TSGetStepNumber(ts, &steps);
81: /* free the memories */
83: TSDestroy(&ts);
84: VecDestroy(&global);
85: MatDestroy(&A);
86: MatDestroy(&S);
88: PetscFinalize();
89: return 0;
90: }
92: PetscErrorCode MyMatMult(Mat S, Vec x, Vec y)
93: {
94: const PetscScalar *inptr;
95: PetscScalar *outptr;
98: VecGetArrayRead(x, &inptr);
99: VecGetArrayWrite(y, &outptr);
101: outptr[0] = 2.0 * inptr[0] + inptr[1];
102: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
103: outptr[2] = inptr[1] + 2.0 * inptr[2];
105: VecRestoreArrayRead(x, &inptr);
106: VecRestoreArrayWrite(y, &outptr);
107: return 0;
108: }
110: /* this test problem has initial values (1,1,1). */
111: PetscErrorCode Initial(Vec global, void *ctx)
112: {
113: PetscScalar *localptr;
114: PetscInt i, mybase, myend, locsize;
116: /* determine starting point of each processor */
117: VecGetOwnershipRange(global, &mybase, &myend);
118: VecGetLocalSize(global, &locsize);
120: /* Initialize the array */
121: VecGetArrayWrite(global, &localptr);
122: for (i = 0; i < locsize; i++) localptr[i] = 1.0;
124: if (mybase == 0) localptr[0] = 1.0;
126: VecRestoreArrayWrite(global, &localptr);
127: return 0;
128: }
130: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx)
131: {
132: VecScatter scatter;
133: IS from, to;
134: PetscInt i, n, *idx;
135: Vec tmp_vec;
136: const PetscScalar *tmp;
138: /* Get the size of the vector */
139: VecGetSize(global, &n);
141: /* Set the index sets */
142: PetscMalloc1(n, &idx);
143: for (i = 0; i < n; i++) idx[i] = i;
145: /* Create local sequential vectors */
146: VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec);
148: /* Create scatter context */
149: ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from);
150: ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to);
151: VecScatterCreate(global, from, tmp_vec, to, &scatter);
152: VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD);
153: VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD);
155: VecGetArrayRead(tmp_vec, &tmp);
156: 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]));
157: 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)));
158: VecRestoreArrayRead(tmp_vec, &tmp);
159: VecScatterDestroy(&scatter);
160: ISDestroy(&from);
161: ISDestroy(&to);
162: PetscFree(idx);
163: VecDestroy(&tmp_vec);
164: return 0;
165: }
167: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
168: {
169: PetscScalar *outptr;
170: const PetscScalar *inptr;
171: PetscInt i, n, *idx;
172: IS from, to;
173: VecScatter scatter;
174: Vec tmp_in, tmp_out;
176: /* Get the length of parallel vector */
177: VecGetSize(globalin, &n);
179: /* Set the index sets */
180: PetscMalloc1(n, &idx);
181: for (i = 0; i < n; i++) idx[i] = i;
183: /* Create local sequential vectors */
184: VecCreateSeq(PETSC_COMM_SELF, n, &tmp_in);
185: VecDuplicate(tmp_in, &tmp_out);
187: /* Create scatter context */
188: ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from);
189: ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to);
190: VecScatterCreate(globalin, from, tmp_in, to, &scatter);
191: VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD);
192: VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD);
193: VecScatterDestroy(&scatter);
195: /*Extract income array */
196: VecGetArrayRead(tmp_in, &inptr);
198: /* Extract outcome array*/
199: VecGetArrayWrite(tmp_out, &outptr);
201: outptr[0] = 2.0 * inptr[0] + inptr[1];
202: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
203: outptr[2] = inptr[1] + 2.0 * inptr[2];
205: VecRestoreArrayRead(tmp_in, &inptr);
206: VecRestoreArrayWrite(tmp_out, &outptr);
208: VecScatterCreate(tmp_out, from, globalout, to, &scatter);
209: VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD);
210: VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD);
212: /* Destroy idx aand scatter */
213: ISDestroy(&from);
214: ISDestroy(&to);
215: VecScatterDestroy(&scatter);
216: VecDestroy(&tmp_in);
217: VecDestroy(&tmp_out);
218: PetscFree(idx);
219: return 0;
220: }
222: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx)
223: {
224: PetscScalar v[3];
225: const PetscScalar *tmp;
226: PetscInt idx[3], i;
228: idx[0] = 0;
229: idx[1] = 1;
230: idx[2] = 2;
231: VecGetArrayRead(x, &tmp);
233: i = 0;
234: v[0] = 2.0;
235: v[1] = 1.0;
236: v[2] = 0.0;
237: MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES);
239: i = 1;
240: v[0] = 1.0;
241: v[1] = 2.0;
242: v[2] = 1.0;
243: MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES);
245: i = 2;
246: v[0] = 0.0;
247: v[1] = 1.0;
248: v[2] = 2.0;
249: MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES);
251: MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY);
252: MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY);
254: if (A != BB) {
255: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
256: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
257: }
258: VecRestoreArrayRead(x, &tmp);
260: return 0;
261: }
263: /*
264: The exact solutions
265: */
266: PetscReal solx(PetscReal t)
267: {
268: 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));
269: }
271: PetscReal soly(PetscReal t)
272: {
273: 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);
274: }
276: PetscReal solz(PetscReal t)
277: {
278: 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));
279: }
281: /*TEST
283: test:
284: suffix: euler
285: args: -ts_type euler
286: requires: !single
288: test:
289: suffix: beuler
290: args: -ts_type beuler
291: requires: !single
293: TEST*/