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