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