1 static char help[] = "Solves DAE with integrator only on non-algebraic terms \n";
2
3 #include <petscts.h>
4
5 /*
6 \dot{U} = f(U,V)
7 F(U,V) = 0
8
9 Same as ex6.c and ex7.c except calls the ARKIMEX integrator on the entire DAE
10 */
11
12 /*
13 f(U,V) = U + V
14
15 */
f(PetscReal t,Vec U,Vec V,Vec F)16 PetscErrorCode f(PetscReal t, Vec U, Vec V, Vec F)
17 {
18 PetscFunctionBeginUser;
19 PetscCall(VecWAXPY(F, 1.0, U, V));
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
23 /*
24 F(U,V) = U - V
25
26 */
F(PetscReal t,Vec U,Vec V,Vec F)27 PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F)
28 {
29 PetscFunctionBeginUser;
30 PetscCall(VecWAXPY(F, -1.0, V, U));
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
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;
41
42 extern PetscErrorCode TSFunctionRHS(TS, PetscReal, Vec, Vec, void *);
43 extern PetscErrorCode TSFunctionI(TS, PetscReal, Vec, Vec, Vec, void *);
44
main(int argc,char ** argv)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;
53
54 PetscFunctionBeginUser;
55 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
56 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
57 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
58 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
59 PetscCall(TSSetType(ts, TSROSW));
60 PetscCall(TSSetFromOptions(ts));
61 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &tsrhs));
62 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &UV));
63 PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunctionRHS, &ctx));
64 PetscCall(TSSetIFunction(ts, NULL, TSFunctionI, &ctx));
65 ctx.f = f;
66 ctx.F = F;
67
68 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.U));
69 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.V));
70 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.UF));
71 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.VF));
72 I = 2 * rank;
73 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is));
74 PetscCall(VecScatterCreate(ctx.U, NULL, UV, is, &ctx.scatterU));
75 PetscCall(ISDestroy(&is));
76 I = 2 * rank + 1;
77 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is));
78 PetscCall(VecScatterCreate(ctx.V, NULL, UV, is, &ctx.scatterV));
79 PetscCall(ISDestroy(&is));
80
81 PetscCall(VecSet(UV, 1.0));
82 PetscCall(TSSolve(ts, UV));
83 PetscCall(VecDestroy(&tsrhs));
84 PetscCall(VecDestroy(&UV));
85 PetscCall(VecDestroy(&ctx.U));
86 PetscCall(VecDestroy(&ctx.V));
87 PetscCall(VecDestroy(&ctx.UF));
88 PetscCall(VecDestroy(&ctx.VF));
89 PetscCall(VecScatterDestroy(&ctx.scatterU));
90 PetscCall(VecScatterDestroy(&ctx.scatterV));
91 PetscCall(TSDestroy(&ts));
92 PetscCall(PetscFinalize());
93 return 0;
94 }
95
96 /*
97 Defines the RHS function that is passed to the time-integrator.
98
99 */
TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void * actx)100 PetscErrorCode TSFunctionRHS(TS ts, PetscReal t, Vec UV, Vec F, void *actx)
101 {
102 AppCtx *ctx = (AppCtx *)actx;
103
104 PetscFunctionBeginUser;
105 PetscCall(VecSet(F, 0.0));
106 PetscCall(VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE));
107 PetscCall(VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE));
108 PetscCall(VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE));
109 PetscCall(VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE));
110 PetscCall((*ctx->f)(t, ctx->U, ctx->V, ctx->UF));
111 PetscCall(VecScatterBegin(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD));
112 PetscCall(VecScatterEnd(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD));
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
116 /*
117 Defines the nonlinear function that is passed to the time-integrator
118
119 */
TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void * actx)120 PetscErrorCode TSFunctionI(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx)
121 {
122 AppCtx *ctx = (AppCtx *)actx;
123
124 PetscFunctionBeginUser;
125 PetscCall(VecCopy(UVdot, F));
126 PetscCall(VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE));
127 PetscCall(VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE));
128 PetscCall(VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE));
129 PetscCall(VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE));
130 PetscCall((*ctx->F)(t, ctx->U, ctx->V, ctx->VF));
131 PetscCall(VecScatterBegin(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD));
132 PetscCall(VecScatterEnd(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD));
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135