Actual source code: ex3fastslowsplit.c
1: static char help[] = "A fast-slow system for testing ARKIMEX.\n";
3: /*F
5: \begin{eqnarray}
6: ys' = -2.0*\frac{-1.0+ys^2.0-\cos(t)}{2.0*ys}+0.05*\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-\frac{\sin(t)}{2.0*ys}\\
7: yf' = 0.05*\frac{-1.0+ys^2-\cos(t)}{2.0*ys}-\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-5.0*\frac{\sin(5.0*t)}{2.0*yf}\\
8: \end{eqnarray}
10: F*/
12: /*
13: This example demonstrates how to use ARKIMEX for solving a fast-slow system. The system is partitioned additively and component-wise at the same time.
14: ys stands for the slow component and yf stands for the fast component. On the RHS for yf, only the term -\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf} is treated implicitly while the rest is treated explicitly.
15: */
17: #include <petscts.h>
19: typedef struct {
20: PetscReal Tf, dt;
21: } AppCtx;
23: static PetscErrorCode RHSFunctionslow(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
24: {
25: const PetscScalar *u;
26: PetscScalar *f;
28: PetscFunctionBegin;
29: PetscCall(VecGetArrayRead(U, &u));
30: PetscCall(VecGetArray(F, &f));
31: f[0] = -2.0 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) + 0.05 * (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - PetscSinScalar(t) / (2.0 * u[0]);
32: PetscCall(VecRestoreArrayRead(U, &u));
33: PetscCall(VecRestoreArray(F, &f));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode RHSFunctionfast(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
38: {
39: const PetscScalar *u;
40: PetscScalar *f;
42: PetscFunctionBegin;
43: PetscCall(VecGetArrayRead(U, &u));
44: PetscCall(VecGetArray(F, &f));
45: f[0] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]);
46: PetscCall(VecRestoreArrayRead(U, &u));
47: PetscCall(VecRestoreArray(F, &f));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode IFunctionfast(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
52: {
53: const PetscScalar *u, *udot;
54: PetscScalar *f;
56: PetscFunctionBegin;
57: PetscCall(VecGetArrayRead(U, &u));
58: PetscCall(VecGetArrayRead(Udot, &udot));
59: PetscCall(VecGetArray(F, &f));
60: f[0] = udot[0] + (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]);
61: PetscCall(VecRestoreArrayRead(Udot, &udot));
62: PetscCall(VecRestoreArrayRead(U, &u));
63: PetscCall(VecRestoreArray(F, &f));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode IJacobianfast(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
68: {
69: PetscInt rowcol[] = {0};
70: const PetscScalar *u;
71: PetscScalar J[1][1];
73: PetscFunctionBeginUser;
74: PetscCall(VecGetArrayRead(U, &u));
75: J[0][0] = a + (2.0 + PetscCosScalar(5.0 * t)) / (2.0 * u[1] * u[1]) + 0.5;
76: PetscCall(MatSetValues(B, 1, rowcol, 1, rowcol, &J[0][0], INSERT_VALUES));
77: PetscCall(VecRestoreArrayRead(U, &u));
79: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
80: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
81: if (A != B) {
82: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
84: }
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: /*
89: Define the analytic solution for check method easily
90: */
91: static PetscErrorCode sol_true(PetscReal t, Vec U)
92: {
93: PetscScalar *u;
95: PetscFunctionBegin;
96: PetscCall(VecGetArray(U, &u));
97: u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t));
98: u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t));
99: PetscCall(VecRestoreArray(U, &u));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: int main(int argc, char **argv)
104: {
105: TS ts; /* ODE integrator */
106: Vec U; /* solution will be stored here */
107: Vec Utrue;
108: Mat A;
109: PetscMPIInt size;
110: AppCtx ctx;
111: PetscScalar *u;
112: IS iss;
113: IS isf;
114: PetscInt *indicess;
115: PetscInt *indicesf;
116: PetscInt n = 2;
117: PetscScalar error, tt;
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Initialize program
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: PetscFunctionBeginUser;
123: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
124: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
125: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Create index for slow part and fast part
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: PetscCall(PetscMalloc1(1, &indicess));
131: indicess[0] = 0;
132: PetscCall(PetscMalloc1(1, &indicesf));
133: indicesf[0] = 1;
134: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss));
135: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf));
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Create necessary vector
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
141: PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
142: PetscCall(VecSetFromOptions(U));
143: PetscCall(VecDuplicate(U, &Utrue));
144: PetscCall(VecCopy(U, Utrue));
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Set runtime options
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
150: ctx.Tf = 0.3;
151: ctx.dt = 0.01;
152: PetscCall(PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL));
153: PetscCall(PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL));
154: PetscOptionsEnd();
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Initialize the solution
158: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: PetscCall(VecGetArray(U, &u));
160: u[0] = PetscSqrtScalar(2.0);
161: u[1] = PetscSqrtScalar(3.0);
162: PetscCall(VecRestoreArray(U, &u));
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Create timestepping solver context
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
168: PetscCall(TSSetType(ts, TSARKIMEX));
169: PetscCall(TSARKIMEXSetFastSlowSplit(ts, PETSC_TRUE));
171: PetscCall(TSRHSSplitSetIS(ts, "slow", iss));
172: PetscCall(TSRHSSplitSetIS(ts, "fast", isf));
173: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunctionFn *)RHSFunctionslow, &ctx));
174: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunctionFn *)RHSFunctionfast, &ctx));
175: PetscCall(TSRHSSplitSetIFunction(ts, "fast", NULL, (TSIFunctionFn *)IFunctionfast, &ctx));
176: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
177: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
178: PetscCall(MatSetFromOptions(A));
179: PetscCall(MatSetUp(A));
180: PetscCall(TSRHSSplitSetIJacobian(ts, "fast", A, A, (TSIJacobianFn *)IJacobianfast, &ctx));
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Set initial conditions
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: PetscCall(TSSetSolution(ts, U));
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Set solver options
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190: PetscCall(TSSetMaxTime(ts, ctx.Tf));
191: PetscCall(TSSetTimeStep(ts, ctx.dt));
192: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
193: PetscCall(TSSetFromOptions(ts));
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Solve linear system
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: PetscCall(TSSolve(ts, U));
199: PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Check the error of the Petsc solution
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: PetscCall(TSGetTime(ts, &tt));
205: PetscCall(sol_true(tt, Utrue));
206: PetscCall(VecAXPY(Utrue, -1.0, U));
207: PetscCall(VecNorm(Utrue, NORM_2, &error));
209: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210: Print norm2 error
211: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)PetscAbsScalar(error)));
214: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: Free work space. All PETSc objects should be destroyed when they are no longer needed.
216: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217: PetscCall(VecDestroy(&U));
218: PetscCall(TSDestroy(&ts));
219: PetscCall(VecDestroy(&Utrue));
220: PetscCall(ISDestroy(&iss));
221: PetscCall(ISDestroy(&isf));
222: PetscCall(PetscFree(indicess));
223: PetscCall(PetscFree(indicesf));
224: PetscCall(MatDestroy(&A));
225: PetscCall(PetscFinalize());
226: return 0;
227: }
229: /*TEST
230: build:
231: requires: !complex
233: test:
235: test:
236: suffix: 2
237: args: -ts_arkimex_initial_guess_extrapolate 1
238: output_file: output/ex3fastslowsplit_1.out
240: TEST*/