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
10 /*
11 f(U,V) = U + V
12 */
f(PetscReal t,Vec U,Vec V,Vec F)13 PetscErrorCode f(PetscReal t, Vec U, Vec V, Vec F)
14 {
15 PetscFunctionBeginUser;
16 PetscCall(VecWAXPY(F, 1.0, U, V));
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
20 /*
21 F(U,V) = U - V
22 */
F(PetscReal t,Vec U,Vec V,Vec F)23 PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F)
24 {
25 PetscFunctionBeginUser;
26 PetscCall(VecWAXPY(F, -1.0, V, U));
27 PetscFunctionReturn(PETSC_SUCCESS);
28 }
29
30 typedef struct {
31 PetscReal t;
32 SNES snes;
33 Vec U, V;
34 PetscErrorCode (*f)(PetscReal, Vec, Vec, Vec);
35 PetscErrorCode (*F)(PetscReal, Vec, Vec, Vec);
36 } AppCtx;
37
38 extern PetscErrorCode TSFunction(TS, PetscReal, Vec, Vec, void *);
39 extern PetscErrorCode SNESFunction(SNES, Vec, Vec, void *);
40
main(int argc,char ** argv)41 int main(int argc, char **argv)
42 {
43 AppCtx ctx;
44 TS ts;
45 Vec tsrhs, U;
46
47 PetscFunctionBeginUser;
48 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
49 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
50 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
51 PetscCall(TSSetType(ts, TSEULER));
52 PetscCall(TSSetFromOptions(ts));
53 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &tsrhs));
54 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &U));
55 PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunction, &ctx));
56 PetscCall(TSSetMaxTime(ts, 1.0));
57 ctx.f = f;
58
59 PetscCall(SNESCreate(PETSC_COMM_WORLD, &ctx.snes));
60 PetscCall(SNESSetFromOptions(ctx.snes));
61 PetscCall(SNESSetFunction(ctx.snes, NULL, SNESFunction, &ctx));
62 PetscCall(SNESSetJacobian(ctx.snes, NULL, NULL, SNESComputeJacobianDefault, &ctx));
63 ctx.F = F;
64 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &ctx.V));
65
66 PetscCall(VecSet(U, 1.0));
67 PetscCall(TSSolve(ts, U));
68
69 PetscCall(VecDestroy(&ctx.V));
70 PetscCall(VecDestroy(&tsrhs));
71 PetscCall(VecDestroy(&U));
72 PetscCall(SNESDestroy(&ctx.snes));
73 PetscCall(TSDestroy(&ts));
74 PetscCall(PetscFinalize());
75 return 0;
76 }
77
78 /*
79 Defines the RHS function that is passed to the time-integrator.
80
81 Solves F(U,V) for V and then computes f(U,V)
82 */
TSFunction(TS ts,PetscReal t,Vec U,Vec F,void * actx)83 PetscErrorCode TSFunction(TS ts, PetscReal t, Vec U, Vec F, void *actx)
84 {
85 AppCtx *ctx = (AppCtx *)actx;
86
87 PetscFunctionBeginUser;
88 ctx->t = t;
89 ctx->U = U;
90 PetscCall(SNESSolve(ctx->snes, NULL, ctx->V));
91 PetscCall((*ctx->f)(t, U, ctx->V, F));
92 PetscFunctionReturn(PETSC_SUCCESS);
93 }
94
95 /*
96 Defines the nonlinear function that is passed to the nonlinear solver
97 */
SNESFunction(SNES snes,Vec V,Vec F,void * actx)98 PetscErrorCode SNESFunction(SNES snes, Vec V, Vec F, void *actx)
99 {
100 AppCtx *ctx = (AppCtx *)actx;
101
102 PetscFunctionBeginUser;
103 PetscCall((*ctx->F)(ctx->t, ctx->U, V, F));
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106
107 /*TEST
108
109 test:
110 args: -ts_monitor -ts_view
111
112 TEST*/
113