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