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;
47: PetscInitialize(&argc,&argv,(char*)0,help);
48: TSCreate(PETSC_COMM_WORLD,&ts);
49: TSSetProblemType(ts,TS_NONLINEAR);
50: TSSetType(ts,TSEULER);
51: TSSetFromOptions(ts);
52: VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs);
53: VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U);
54: TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
55: TSSetMaxTime(ts,1.0);
56: ctx.f = f;
58: SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
59: SNESSetFromOptions(ctx.snes);
60: SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
61: SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
62: ctx.F = F;
63: VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V);
65: VecSet(U,1.0);
66: TSSolve(ts,U);
68: VecDestroy(&ctx.V);
69: VecDestroy(&tsrhs);
70: VecDestroy(&U);
71: SNESDestroy(&ctx.snes);
72: TSDestroy(&ts);
73: PetscFinalize();
74: return 0;
75: }
77: /*
78: Defines the RHS function that is passed to the time-integrator.
80: Solves F(U,V) for V and then computes f(U,V)
81: */
82: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
83: {
84: AppCtx *ctx = (AppCtx*)actx;
87: ctx->t = t;
88: ctx->U = U;
89: SNESSolve(ctx->snes,NULL,ctx->V);
90: (*ctx->f)(t,U,ctx->V,F);
91: return 0;
92: }
94: /*
95: Defines the nonlinear function that is passed to the nonlinear solver
96: */
97: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
98: {
99: AppCtx *ctx = (AppCtx*)actx;
102: (*ctx->F)(ctx->t,ctx->U,V,F);
103: return 0;
104: }
106: /*TEST
108: test:
109: args: -ts_monitor -ts_view
111: TEST*/