1 static char help[] = "Solves DAE with integrator only on non-algebraic terms \n";
2
3 #include <petscts.h>
4
5 /*
6 \dot{U} = f(U,V)
7 F(U,V) = 0
8 */
9
10
11 /*
12 f(U,V) = U + V
13 */
f(PetscReal t,Vec U,Vec V,Vec F)14 PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F)
15 {
16 PetscErrorCode ierr;
17
18 PetscFunctionBeginUser;
19 ierr = VecWAXPY(F,1.0,U,V);CHKERRQ(ierr);
20 PetscFunctionReturn(0);
21 }
22
23 /*
24 F(U,V) = U - V
25 */
F(PetscReal t,Vec U,Vec V,Vec F)26 PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
27 {
28 PetscErrorCode ierr;
29
30 PetscFunctionBeginUser;
31 ierr = VecWAXPY(F,-1.0,V,U);CHKERRQ(ierr);
32 PetscFunctionReturn(0);
33 }
34
35 typedef struct {
36 PetscReal t;
37 SNES snes;
38 Vec U,V;
39 PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
40 PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
41 } AppCtx;
42
43 extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
44 extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);
45
main(int argc,char ** argv)46 int main(int argc,char **argv)
47 {
48 PetscErrorCode ierr;
49 AppCtx ctx;
50 TS ts;
51 Vec tsrhs,U;
52
53 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
54 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
55 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
56 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr);
57 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
58 ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs);CHKERRQ(ierr);
59 ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U);CHKERRQ(ierr);
60 ierr = TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);CHKERRQ(ierr);
61 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr);
62 ctx.f = f;
63
64 ierr = SNESCreate(PETSC_COMM_WORLD,&ctx.snes);CHKERRQ(ierr);
65 ierr = SNESSetFromOptions(ctx.snes);CHKERRQ(ierr);
66 ierr = SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);CHKERRQ(ierr);
67 ierr = SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);CHKERRQ(ierr);
68 ctx.F = F;
69 ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V);CHKERRQ(ierr);
70
71 ierr = VecSet(U,1.0);CHKERRQ(ierr);
72 ierr = TSSolve(ts,U);CHKERRQ(ierr);
73
74 ierr = VecDestroy(&ctx.V);CHKERRQ(ierr);
75 ierr = VecDestroy(&tsrhs);CHKERRQ(ierr);
76 ierr = VecDestroy(&U);CHKERRQ(ierr);
77 ierr = SNESDestroy(&ctx.snes);CHKERRQ(ierr);
78 ierr = TSDestroy(&ts);CHKERRQ(ierr);
79 ierr = PetscFinalize();
80 return ierr;
81 }
82
83 /*
84 Defines the RHS function that is passed to the time-integrator.
85
86 Solves F(U,V) for V and then computes f(U,V)
87 */
TSFunction(TS ts,PetscReal t,Vec U,Vec F,void * actx)88 PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
89 {
90 AppCtx *ctx = (AppCtx*)actx;
91 PetscErrorCode ierr;
92
93 PetscFunctionBeginUser;
94 ctx->t = t;
95 ctx->U = U;
96 ierr = SNESSolve(ctx->snes,NULL,ctx->V);CHKERRQ(ierr);
97 ierr = (*ctx->f)(t,U,ctx->V,F);CHKERRQ(ierr);
98 PetscFunctionReturn(0);
99 }
100
101 /*
102 Defines the nonlinear function that is passed to the nonlinear solver
103 */
SNESFunction(SNES snes,Vec V,Vec F,void * actx)104 PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
105 {
106 AppCtx *ctx = (AppCtx*)actx;
107 PetscErrorCode ierr;
108
109 PetscFunctionBeginUser;
110 ierr = (*ctx->F)(ctx->t,ctx->U,V,F);CHKERRQ(ierr);
111 PetscFunctionReturn(0);
112 }
113
114 /*TEST
115
116 test:
117 args: -ts_monitor -ts_view
118
119 TEST*/
120