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