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 TSROSW integrator on the entire DAE
10 */
11
12 /*
13 f(U,V) = U + V
14
15 */
f(PetscReal t,Vec UV,Vec F)16 PetscErrorCode f(PetscReal t, Vec UV, Vec F)
17 {
18 const PetscScalar *u, *v;
19 PetscScalar *f;
20 PetscInt n, i;
21
22 PetscFunctionBeginUser;
23 PetscCall(VecGetLocalSize(UV, &n));
24 n = n / 2;
25 PetscCall(VecGetArrayRead(UV, &u));
26 v = u + n;
27 PetscCall(VecGetArrayWrite(F, &f));
28 for (i = 0; i < n; i++) f[i] = u[i] + v[i];
29 PetscCall(VecRestoreArrayRead(UV, &u));
30 PetscCall(VecRestoreArrayWrite(F, &f));
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
34 /*
35 F(U,V) = U - V
36
37 */
F(PetscReal t,Vec UV,Vec F)38 PetscErrorCode F(PetscReal t, Vec UV, Vec F)
39 {
40 const PetscScalar *u, *v;
41 PetscScalar *f;
42 PetscInt n, i;
43
44 PetscFunctionBeginUser;
45 PetscCall(VecGetLocalSize(UV, &n));
46 n = n / 2;
47 PetscCall(VecGetArrayRead(UV, &u));
48 v = u + n;
49 PetscCall(VecGetArrayWrite(F, &f));
50 f = f + n;
51 for (i = 0; i < n; i++) f[i] = u[i] - v[i];
52 f = f - n;
53 PetscCall(VecRestoreArrayRead(UV, &u));
54 PetscCall(VecRestoreArrayWrite(F, &f));
55 PetscFunctionReturn(PETSC_SUCCESS);
56 }
57
58 typedef struct {
59 PetscErrorCode (*f)(PetscReal, Vec, Vec);
60 PetscErrorCode (*F)(PetscReal, Vec, Vec);
61 } AppCtx;
62
63 extern PetscErrorCode TSFunctionRHS(TS, PetscReal, Vec, Vec, void *);
64 extern PetscErrorCode TSFunctionI(TS, PetscReal, Vec, Vec, Vec, void *);
65
main(int argc,char ** argv)66 int main(int argc, char **argv)
67 {
68 AppCtx ctx;
69 TS ts;
70 Vec tsrhs, UV;
71
72 PetscFunctionBeginUser;
73 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
74 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
75 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
76 PetscCall(TSSetType(ts, TSROSW));
77 PetscCall(TSSetFromOptions(ts));
78 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &tsrhs));
79 PetscCall(VecDuplicate(tsrhs, &UV));
80 PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunctionRHS, &ctx));
81 PetscCall(TSSetIFunction(ts, NULL, TSFunctionI, &ctx));
82 PetscCall(TSSetMaxTime(ts, 1.0));
83 ctx.f = f;
84 ctx.F = F;
85
86 PetscCall(VecSet(UV, 1.0));
87 PetscCall(TSSolve(ts, UV));
88 PetscCall(VecDestroy(&tsrhs));
89 PetscCall(VecDestroy(&UV));
90 PetscCall(TSDestroy(&ts));
91 PetscCall(PetscFinalize());
92 return 0;
93 }
94
95 /*
96 Defines the RHS function that is passed to the time-integrator.
97 */
TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void * actx)98 PetscErrorCode TSFunctionRHS(TS ts, PetscReal t, Vec UV, Vec F, void *actx)
99 {
100 AppCtx *ctx = (AppCtx *)actx;
101
102 PetscFunctionBeginUser;
103 PetscCall(VecSet(F, 0.0));
104 PetscCall((*ctx->f)(t, UV, F));
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
108 /*
109 Defines the nonlinear function that is passed to the time-integrator
110 */
TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void * actx)111 PetscErrorCode TSFunctionI(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx)
112 {
113 AppCtx *ctx = (AppCtx *)actx;
114
115 PetscFunctionBeginUser;
116 PetscCall(VecCopy(UVdot, F));
117 PetscCall((*ctx->F)(t, UV, F));
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
121 /*TEST
122
123 test:
124 args: -ts_view
125
126 test:
127 suffix: 2
128 args: -snes_lag_jacobian 2 -ts_view
129
130 TEST*/
131