xref: /petsc/src/ts/tests/ex6.c (revision 5f80ce2ab25dff0f4601e710601cbbcecf323266) 
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*5f80ce2aSJacob Faibussowitsch   CHKERRQ(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*5f80ce2aSJacob Faibussowitsch   CHKERRQ(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   PetscErrorCode ierr;
44c4762a1bSJed Brown   AppCtx         ctx;
45c4762a1bSJed Brown   TS             ts;
46c4762a1bSJed Brown   Vec            tsrhs,U;
47c4762a1bSJed Brown 
48c4762a1bSJed Brown   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
49*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts));
50*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR));
51*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetType(ts,TSEULER));
52*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetFromOptions(ts));
53*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs));
54*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U));
55*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx));
56*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSetMaxTime(ts,1.0));
57c4762a1bSJed Brown   ctx.f = f;
58c4762a1bSJed Brown 
59*5f80ce2aSJacob Faibussowitsch   CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&ctx.snes));
60*5f80ce2aSJacob Faibussowitsch   CHKERRQ(SNESSetFromOptions(ctx.snes));
61*5f80ce2aSJacob Faibussowitsch   CHKERRQ(SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx));
62*5f80ce2aSJacob Faibussowitsch   CHKERRQ(SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx));
63c4762a1bSJed Brown   ctx.F = F;
64*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V));
65c4762a1bSJed Brown 
66*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(U,1.0));
67*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSSolve(ts,U));
68c4762a1bSJed Brown 
69*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&ctx.V));
70*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&tsrhs));
71*5f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&U));
72*5f80ce2aSJacob Faibussowitsch   CHKERRQ(SNESDestroy(&ctx.snes));
73*5f80ce2aSJacob Faibussowitsch   CHKERRQ(TSDestroy(&ts));
74c4762a1bSJed Brown   ierr = PetscFinalize();
75c4762a1bSJed Brown   return ierr;
76c4762a1bSJed Brown }
77c4762a1bSJed Brown 
78c4762a1bSJed Brown /*
79c4762a1bSJed Brown    Defines the RHS function that is passed to the time-integrator.
80c4762a1bSJed Brown 
81c4762a1bSJed Brown    Solves F(U,V) for V and then computes f(U,V)
82c4762a1bSJed Brown */
83c4762a1bSJed Brown PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
84c4762a1bSJed Brown {
85c4762a1bSJed Brown   AppCtx         *ctx = (AppCtx*)actx;
86c4762a1bSJed Brown 
87c4762a1bSJed Brown   PetscFunctionBeginUser;
88c4762a1bSJed Brown   ctx->t = t;
89c4762a1bSJed Brown   ctx->U = U;
90*5f80ce2aSJacob Faibussowitsch   CHKERRQ(SNESSolve(ctx->snes,NULL,ctx->V));
91*5f80ce2aSJacob Faibussowitsch   CHKERRQ((*ctx->f)(t,U,ctx->V,F));
92c4762a1bSJed Brown   PetscFunctionReturn(0);
93c4762a1bSJed Brown }
94c4762a1bSJed Brown 
95c4762a1bSJed Brown /*
96c4762a1bSJed Brown    Defines the nonlinear function that is passed to the nonlinear solver
97c4762a1bSJed Brown */
98c4762a1bSJed Brown PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
99c4762a1bSJed Brown {
100c4762a1bSJed Brown   AppCtx         *ctx = (AppCtx*)actx;
101c4762a1bSJed Brown 
102c4762a1bSJed Brown   PetscFunctionBeginUser;
103*5f80ce2aSJacob Faibussowitsch   CHKERRQ((*ctx->F)(ctx->t,ctx->U,V,F));
104c4762a1bSJed Brown   PetscFunctionReturn(0);
105c4762a1bSJed Brown }
106c4762a1bSJed Brown 
107303a5415SBarry Smith /*TEST
108c4762a1bSJed Brown 
109303a5415SBarry Smith     test:
110303a5415SBarry Smith       args:  -ts_monitor -ts_view
111303a5415SBarry Smith 
112303a5415SBarry Smith TEST*/
113