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]
11: When run in parallel, each process solves the same set of equations separately.
12: */
14: static char help[] = "Solves a linear ODE. \n\n";
16: #include <petscts.h>
17: #include <petscpc.h>
19: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
20: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
21: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
22: extern PetscErrorCode Initial(Vec, void *);
23: extern PetscErrorCode MyMatMult(Mat, Vec, Vec);
25: extern PetscReal solx(PetscReal);
26: extern PetscReal soly(PetscReal);
27: extern PetscReal solz(PetscReal);
29: int main(int argc, char **argv)
30: {
31: PetscInt time_steps = 100, steps;
32: Vec global;
33: PetscReal dt, ftime;
34: TS ts;
35: Mat A, S;
36: PetscBool nest = PETSC_FALSE;
38: PetscFunctionBeginUser;
39: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
40: PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
41: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nest", &nest, NULL));
43: /* create vector to hold state */
44: if (nest) {
45: Vec g[3];
47: PetscCall(VecCreate(PETSC_COMM_WORLD, &g[0]));
48: PetscCall(VecSetSizes(g[0], 1, PETSC_DECIDE));
49: PetscCall(VecSetFromOptions(g[0]));
50: PetscCall(VecDuplicate(g[0], &g[1]));
51: PetscCall(VecDuplicate(g[0], &g[2]));
52: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 3, NULL, g, &global));
53: PetscCall(VecDestroy(&g[0]));
54: PetscCall(VecDestroy(&g[1]));
55: PetscCall(VecDestroy(&g[2]));
56: } else {
57: PetscCall(VecCreate(PETSC_COMM_WORLD, &global));
58: PetscCall(VecSetSizes(global, 3, PETSC_DECIDE));
59: PetscCall(VecSetFromOptions(global));
60: }
62: /* set initial conditions */
63: PetscCall(Initial(global, NULL));
65: /* make timestep context */
66: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
67: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
68: PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
69: dt = 0.001;
71: /*
72: The user provides the RHS and Jacobian
73: */
74: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
75: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
76: PetscCall(MatSetSizes(A, 3, 3, PETSC_DECIDE, PETSC_DECIDE));
77: PetscCall(MatSetFromOptions(A));
78: PetscCall(MatSetUp(A));
79: PetscCall(RHSJacobian(ts, 0.0, global, A, A, NULL));
80: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
82: PetscCall(MatCreateShell(PETSC_COMM_WORLD, 3, 3, PETSC_DECIDE, PETSC_DECIDE, NULL, &S));
83: PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MyMatMult));
84: PetscCall(TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL));
86: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
87: PetscCall(TSSetFromOptions(ts));
89: PetscCall(TSSetTimeStep(ts, dt));
90: PetscCall(TSSetMaxSteps(ts, time_steps));
91: PetscCall(TSSetMaxTime(ts, 1));
92: PetscCall(TSSetSolution(ts, global));
94: PetscCall(TSSolve(ts, global));
95: PetscCall(TSGetSolveTime(ts, &ftime));
96: PetscCall(TSGetStepNumber(ts, &steps));
98: /* free the memory */
99: PetscCall(TSDestroy(&ts));
100: PetscCall(VecDestroy(&global));
101: PetscCall(MatDestroy(&A));
102: PetscCall(MatDestroy(&S));
104: PetscCall(PetscFinalize());
105: return 0;
106: }
108: PetscErrorCode MyMatMult(Mat S, Vec x, Vec y)
109: {
110: const PetscScalar *inptr;
111: PetscScalar *outptr;
113: PetscFunctionBeginUser;
114: PetscCall(VecGetArrayRead(x, &inptr));
115: PetscCall(VecGetArrayWrite(y, &outptr));
117: outptr[0] = 2.0 * inptr[0] + inptr[1];
118: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
119: outptr[2] = inptr[1] + 2.0 * inptr[2];
121: PetscCall(VecRestoreArrayRead(x, &inptr));
122: PetscCall(VecRestoreArrayWrite(y, &outptr));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: PetscErrorCode Initial(Vec global, void *ctx)
127: {
128: PetscScalar *localptr;
130: PetscFunctionBeginUser;
131: PetscCall(VecGetArrayWrite(global, &localptr));
132: localptr[0] = solx(0.0);
133: localptr[1] = soly(0.0);
134: localptr[2] = solz(0.0);
135: PetscCall(VecRestoreArrayWrite(global, &localptr));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx)
140: {
141: const PetscScalar *tmp;
142: PetscScalar exact[] = {solx(time), soly(time), solz(time)};
144: PetscFunctionBeginUser;
145: PetscCall(VecGetArrayRead(global, &tmp));
146: 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])));
147: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e errors = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0] - exact[0]), (double)PetscRealPart(tmp[1] - exact[1]), (double)PetscRealPart(tmp[2] - exact[2])));
148: PetscCall(VecRestoreArrayRead(global, &tmp));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
153: {
154: PetscScalar *outptr;
155: const PetscScalar *inptr;
157: PetscFunctionBeginUser;
158: /*Extract income array */
159: PetscCall(VecGetArrayRead(globalin, &inptr));
161: /* Extract outcome array*/
162: PetscCall(VecGetArrayWrite(globalout, &outptr));
164: outptr[0] = 2.0 * inptr[0] + inptr[1];
165: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
166: outptr[2] = inptr[1] + 2.0 * inptr[2];
168: PetscCall(VecRestoreArrayRead(globalin, &inptr));
169: PetscCall(VecRestoreArrayWrite(globalout, &outptr));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx)
174: {
175: PetscScalar v[3];
176: PetscInt idx[3], rst;
178: PetscFunctionBeginUser;
179: PetscCall(VecGetOwnershipRange(x, &rst, NULL));
180: idx[0] = 0 + rst;
181: idx[1] = 1 + rst;
182: idx[2] = 2 + rst;
184: v[0] = 2.0;
185: v[1] = 1.0;
186: v[2] = 0.0;
187: PetscCall(MatSetValues(BB, 1, idx, 3, idx, v, INSERT_VALUES));
189: v[0] = 1.0;
190: v[1] = 2.0;
191: v[2] = 1.0;
192: PetscCall(MatSetValues(BB, 1, idx + 1, 3, idx, v, INSERT_VALUES));
194: v[0] = 0.0;
195: v[1] = 1.0;
196: v[2] = 2.0;
197: PetscCall(MatSetValues(BB, 1, idx + 2, 3, idx, v, INSERT_VALUES));
199: PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
200: PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));
202: if (A != BB) {
203: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
204: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
205: }
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*
210: The exact solutions
211: */
212: PetscReal solx(PetscReal t)
213: {
214: 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));
215: }
217: PetscReal soly(PetscReal t)
218: {
219: 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);
220: }
222: PetscReal solz(PetscReal t)
223: {
224: 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));
225: }
227: /*TEST
229: test:
230: suffix: euler
231: args: -ts_type euler -nest {{0 1}}
232: requires: !single
234: test:
235: suffix: beuler
236: args: -ts_type beuler -nest {{0 1}}
237: requires: !single
239: test:
240: suffix: rk
241: args: -ts_type rk -nest {{0 1}} -ts_adapt_monitor
242: requires: !single
244: test:
245: diff_args: -j
246: requires: double !complex
247: output_file: output/ex2_be_adapt.out
248: suffix: bdf_1_adapt
249: args: -ts_type bdf -ts_bdf_order 1 -ts_adapt_type basic -ts_adapt_clip 0,2
251: test:
252: diff_args: -j
253: requires: double !complex
254: suffix: be_adapt
255: args: -ts_type beuler -ts_adapt_type basic -ts_adapt_clip 0,2
257: TEST*/