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;
72: PetscInitialize(&argc,&argv,(char*)0,help);
73: TSCreate(PETSC_COMM_WORLD,&ts);
74: TSSetProblemType(ts,TS_NONLINEAR);
75: TSSetType(ts,TSROSW);
76: TSSetFromOptions(ts);
77: VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);
78: VecDuplicate(tsrhs,&UV);
79: TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);
80: TSSetIFunction(ts,NULL,TSFunctionI,&ctx);
81: TSSetMaxTime(ts,1.0);
82: ctx.f = f;
83: ctx.F = F;
85: VecSet(UV,1.0);
86: TSSolve(ts,UV);
87: VecDestroy(&tsrhs);
88: VecDestroy(&UV);
89: TSDestroy(&ts);
90: PetscFinalize();
91: return 0;
92: }
94: /*
95: Defines the RHS function that is passed to the time-integrator.
96: */
97: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
98: {
99: AppCtx *ctx = (AppCtx*)actx;
102: VecSet(F,0.0);
103: (*ctx->f)(t,UV,F);
104: return 0;
105: }
107: /*
108: Defines the nonlinear function that is passed to the time-integrator
109: */
110: PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
111: {
112: AppCtx *ctx = (AppCtx*)actx;
115: VecCopy(UVdot,F);
116: (*ctx->F)(t,UV,F);
117: return 0;
118: }
120: /*TEST
122: test:
123: args: -ts_view
125: test:
126: suffix: 2
127: args: -snes_lag_jacobian 2 -ts_view
129: TEST*/