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;
37: TSType tstype;
38: const char **types;
39: int ntypes;
40: char ststype[16];
42: PetscFunctionBeginUser;
43: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
44: PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
45: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nest", &nest, NULL));
47: /* create vector to hold state */
48: if (nest) {
49: Vec g[3];
51: PetscCall(VecCreate(PETSC_COMM_WORLD, &g[0]));
52: PetscCall(VecSetSizes(g[0], 1, PETSC_DECIDE));
53: PetscCall(VecSetFromOptions(g[0]));
54: PetscCall(VecDuplicate(g[0], &g[1]));
55: PetscCall(VecDuplicate(g[0], &g[2]));
56: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 3, NULL, g, &global));
57: PetscCall(VecDestroy(&g[0]));
58: PetscCall(VecDestroy(&g[1]));
59: PetscCall(VecDestroy(&g[2]));
60: } else {
61: PetscCall(VecCreate(PETSC_COMM_WORLD, &global));
62: PetscCall(VecSetSizes(global, 3, PETSC_DECIDE));
63: PetscCall(VecSetFromOptions(global));
64: }
66: /* set initial conditions */
67: PetscCall(Initial(global, NULL));
69: /* make timestep context */
70: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
71: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
72: PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
73: dt = 0.001;
75: /*
76: The user provides the RHS and Jacobian
77: */
78: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
79: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
80: PetscCall(MatSetSizes(A, 3, 3, PETSC_DECIDE, PETSC_DECIDE));
81: PetscCall(MatSetFromOptions(A));
82: PetscCall(MatSetUp(A));
83: PetscCall(RHSJacobian(ts, 0.0, global, A, A, NULL));
84: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
86: PetscCall(MatCreateShell(PETSC_COMM_WORLD, 3, 3, PETSC_DECIDE, PETSC_DECIDE, NULL, &S));
87: PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MyMatMult));
88: PetscCall(TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL));
90: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
91: PetscCall(TSSetFromOptions(ts));
93: PetscCall(TSSetTimeStep(ts, dt));
94: PetscCall(TSSetMaxSteps(ts, time_steps));
95: PetscCall(TSSetMaxTime(ts, 1));
96: PetscCall(TSSetSolution(ts, global));
98: PetscCall(TSSetUp(ts));
99: PetscCall(TSGetType(ts, &tstype));
100: PetscCall(PetscStrncpy(ststype, tstype, sizeof(ststype)));
101: PetscCall(PetscFunctionListGet(TSList, &types, &ntypes));
102: for (PetscInt i = 0; i < ntypes; i++) {
103: PetscBool reject;
105: PetscCall(TSSetType(ts, types[i]));
106: PetscCall(PetscObjectTypeCompareAny((PetscObject)ts, &reject, TSIRK, TSBASICSYMPLECTIC, TSMPRK, ""));
107: if (!reject) PetscCall(TSSetUp(ts));
108: }
109: PetscCall(PetscFree(types));
110: PetscCall(TSSetType(ts, ststype));
111: PetscCall(TSSetFromOptions(ts));
113: PetscCall(TSSolve(ts, global));
114: PetscCall(TSGetSolveTime(ts, &ftime));
115: PetscCall(TSGetStepNumber(ts, &steps));
117: /* free the memory */
118: PetscCall(TSDestroy(&ts));
119: PetscCall(VecDestroy(&global));
120: PetscCall(MatDestroy(&A));
121: PetscCall(MatDestroy(&S));
123: PetscCall(PetscFinalize());
124: return 0;
125: }
127: PetscErrorCode MyMatMult(Mat S, Vec x, Vec y)
128: {
129: const PetscScalar *inptr;
130: PetscScalar *outptr;
132: PetscFunctionBeginUser;
133: PetscCall(VecGetArrayRead(x, &inptr));
134: PetscCall(VecGetArrayWrite(y, &outptr));
136: outptr[0] = 2.0 * inptr[0] + inptr[1];
137: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
138: outptr[2] = inptr[1] + 2.0 * inptr[2];
140: PetscCall(VecRestoreArrayRead(x, &inptr));
141: PetscCall(VecRestoreArrayWrite(y, &outptr));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: PetscErrorCode Initial(Vec global, PetscCtx ctx)
146: {
147: PetscScalar *localptr;
149: PetscFunctionBeginUser;
150: PetscCall(VecGetArrayWrite(global, &localptr));
151: localptr[0] = solx(0.0);
152: localptr[1] = soly(0.0);
153: localptr[2] = solz(0.0);
154: PetscCall(VecRestoreArrayWrite(global, &localptr));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, PetscCtx ctx)
159: {
160: const PetscScalar *tmp;
161: PetscScalar exact[] = {solx(time), soly(time), solz(time)};
163: PetscFunctionBeginUser;
164: PetscCall(VecGetArrayRead(global, &tmp));
165: 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])));
166: 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])));
167: PetscCall(VecRestoreArrayRead(global, &tmp));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, PetscCtx ctx)
172: {
173: PetscScalar *outptr;
174: const PetscScalar *inptr;
176: PetscFunctionBeginUser;
177: /*Extract income array */
178: PetscCall(VecGetArrayRead(globalin, &inptr));
180: /* Extract outcome array*/
181: PetscCall(VecGetArrayWrite(globalout, &outptr));
183: outptr[0] = 2.0 * inptr[0] + inptr[1];
184: outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
185: outptr[2] = inptr[1] + 2.0 * inptr[2];
187: PetscCall(VecRestoreArrayRead(globalin, &inptr));
188: PetscCall(VecRestoreArrayWrite(globalout, &outptr));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, PetscCtx ctx)
193: {
194: PetscScalar v[3];
195: PetscInt idx[3], rst;
197: PetscFunctionBeginUser;
198: PetscCall(VecGetOwnershipRange(x, &rst, NULL));
199: idx[0] = 0 + rst;
200: idx[1] = 1 + rst;
201: idx[2] = 2 + rst;
203: v[0] = 2.0;
204: v[1] = 1.0;
205: v[2] = 0.0;
206: PetscCall(MatSetValues(BB, 1, idx, 3, idx, v, INSERT_VALUES));
208: v[0] = 1.0;
209: v[1] = 2.0;
210: v[2] = 1.0;
211: PetscCall(MatSetValues(BB, 1, idx + 1, 3, idx, v, INSERT_VALUES));
213: v[0] = 0.0;
214: v[1] = 1.0;
215: v[2] = 2.0;
216: PetscCall(MatSetValues(BB, 1, idx + 2, 3, idx, v, INSERT_VALUES));
218: PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
219: PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));
221: if (A != BB) {
222: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
223: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
224: }
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*
229: The exact solutions
230: */
231: PetscReal solx(PetscReal t)
232: {
233: 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));
234: }
236: PetscReal soly(PetscReal t)
237: {
238: 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);
239: }
241: PetscReal solz(PetscReal t)
242: {
243: 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));
244: }
246: /*TEST
248: test:
249: suffix: euler
250: args: -ts_type euler -nest {{0 1}}
251: requires: !single !sundials2
253: test:
254: suffix: beuler
255: args: -ts_type beuler -nest {{0 1}}
256: requires: !single !sundials2
258: test:
259: suffix: rk
260: args: -ts_type rk -nest {{0 1}} -ts_adapt_monitor
261: requires: !single !sundials2
263: test:
264: diff_args: -j
265: requires: double !complex !sundials2
266: output_file: output/ex2_be_adapt.out
267: suffix: bdf_1_adapt
268: args: -ts_type bdf -ts_bdf_order 1 -ts_adapt_type basic -ts_adapt_clip 0,2
270: test:
271: diff_args: -j
272: requires: double !complex !sundials2
273: suffix: be_adapt
274: args: -ts_type beuler -ts_adapt_type basic -ts_adapt_clip 0,2
276: TEST*/