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