Actual source code: ex9.c
1: static char help[] = "Solves DAE with integrator only on non-algebraic terms \n";
3: #include <petscts.h>
5: /*
6: \dot{U} = f(U,V)
7: F(U,V) = 0
9: Same as ex6.c and ex7.c except calls the ARKIMEX integrator on the entire DAE
10: */
12: /*
13: f(U,V) = U + V
15: */
16: PetscErrorCode f(PetscReal t, Vec U, Vec V, Vec F)
17: {
19: VecWAXPY(F, 1.0, U, V);
20: return 0;
21: }
23: /*
24: F(U,V) = U - V
26: */
27: PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F)
28: {
30: VecWAXPY(F, -1.0, V, U);
31: return 0;
32: }
34: typedef struct {
35: Vec U, V;
36: Vec UF, VF;
37: VecScatter scatterU, scatterV;
38: PetscErrorCode (*f)(PetscReal, Vec, Vec, Vec);
39: PetscErrorCode (*F)(PetscReal, Vec, Vec, Vec);
40: } AppCtx;
42: extern PetscErrorCode TSFunctionRHS(TS, PetscReal, Vec, Vec, void *);
43: extern PetscErrorCode TSFunctionI(TS, PetscReal, Vec, Vec, Vec, void *);
45: int main(int argc, char **argv)
46: {
47: AppCtx ctx;
48: TS ts;
49: Vec tsrhs, UV;
50: IS is;
51: PetscInt I;
52: PetscMPIInt rank;
55: PetscInitialize(&argc, &argv, (char *)0, help);
56: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
57: TSCreate(PETSC_COMM_WORLD, &ts);
58: TSSetProblemType(ts, TS_NONLINEAR);
59: TSSetType(ts, TSROSW);
60: TSSetFromOptions(ts);
61: VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &tsrhs);
62: VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &UV);
63: TSSetRHSFunction(ts, tsrhs, TSFunctionRHS, &ctx);
64: TSSetIFunction(ts, NULL, TSFunctionI, &ctx);
65: ctx.f = f;
66: ctx.F = F;
68: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.U);
69: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.V);
70: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.UF);
71: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.VF);
72: I = 2 * rank;
73: ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is);
74: VecScatterCreate(ctx.U, NULL, UV, is, &ctx.scatterU);
75: ISDestroy(&is);
76: I = 2 * rank + 1;
77: ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is);
78: VecScatterCreate(ctx.V, NULL, UV, is, &ctx.scatterV);
79: ISDestroy(&is);
81: VecSet(UV, 1.0);
82: TSSolve(ts, UV);
83: VecDestroy(&tsrhs);
84: VecDestroy(&UV);
85: VecDestroy(&ctx.U);
86: VecDestroy(&ctx.V);
87: VecDestroy(&ctx.UF);
88: VecDestroy(&ctx.VF);
89: VecScatterDestroy(&ctx.scatterU);
90: VecScatterDestroy(&ctx.scatterV);
91: TSDestroy(&ts);
92: PetscFinalize();
93: return 0;
94: }
96: /*
97: Defines the RHS function that is passed to the time-integrator.
99: */
100: PetscErrorCode TSFunctionRHS(TS ts, PetscReal t, Vec UV, Vec F, void *actx)
101: {
102: AppCtx *ctx = (AppCtx *)actx;
105: VecSet(F, 0.0);
106: VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
107: VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
108: VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
109: VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
110: (*ctx->f)(t, ctx->U, ctx->V, ctx->UF);
111: VecScatterBegin(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD);
112: VecScatterEnd(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD);
113: return 0;
114: }
116: /*
117: Defines the nonlinear function that is passed to the time-integrator
119: */
120: PetscErrorCode TSFunctionI(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx)
121: {
122: AppCtx *ctx = (AppCtx *)actx;
125: VecCopy(UVdot, F);
126: VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
127: VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
128: VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
129: VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
130: (*ctx->F)(t, ctx->U, ctx->V, ctx->VF);
131: VecScatterBegin(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD);
132: VecScatterEnd(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD);
133: return 0;
134: }