Actual source code: ex9.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 ARKIMEX integrator on the entire DAE
10: */
12: /*
13: f(U,V) = U + V
15: */
16: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F)
17: {
19: VecWAXPY(F,1.0,U,V);
20: return 0;
21: }
23: /*
24: F(U,V) = U - V
26: */
27: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F)
28: {
30: VecWAXPY(F,-1.0,V,U);
31: return 0;
32: }
34: typedef struct {
35: Vec U,V;
36: Vec UF,VF;
37: VecScatter scatterU,scatterV;
38: PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec);
39: PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec);
40: } AppCtx;
42: extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*);
43: extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*);
45: int main(int argc,char **argv)
46: {
47: AppCtx ctx;
48: TS ts;
49: Vec tsrhs,UV;
50: IS is;
51: PetscInt I;
52: PetscMPIInt rank;
54: PetscInitialize(&argc,&argv,(char*)0,help);
55: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
56: TSCreate(PETSC_COMM_WORLD,&ts);
57: TSSetProblemType(ts,TS_NONLINEAR);
58: TSSetType(ts,TSROSW);
59: TSSetFromOptions(ts);
60: VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs);
61: VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV);
62: TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);
63: TSSetIFunction(ts,NULL,TSFunctionI,&ctx);
64: ctx.f = f;
65: ctx.F = F;
67: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.U);
68: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);
69: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.UF);
70: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.VF);
71: I = 2*rank;
72: ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
73: VecScatterCreate(ctx.U,NULL,UV,is,&ctx.scatterU);
74: ISDestroy(&is);
75: I = 2*rank + 1;
76: ISCreateGeneral(PETSC_COMM_WORLD,1,&I,PETSC_COPY_VALUES,&is);
77: VecScatterCreate(ctx.V,NULL,UV,is,&ctx.scatterV);
78: ISDestroy(&is);
80: VecSet(UV,1.0);
81: TSSolve(ts,UV);
82: VecDestroy(&tsrhs);
83: VecDestroy(&UV);
84: VecDestroy(&ctx.U);
85: VecDestroy(&ctx.V);
86: VecDestroy(&ctx.UF);
87: VecDestroy(&ctx.VF);
88: VecScatterDestroy(&ctx.scatterU);
89: VecScatterDestroy(&ctx.scatterV);
90: TSDestroy(&ts);
91: PetscFinalize();
92: return 0;
93: }
95: /*
96: Defines the RHS function that is passed to the time-integrator.
98: */
99: PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
100: {
101: AppCtx *ctx = (AppCtx*)actx;
104: VecSet(F,0.0);
105: VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
106: VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
107: VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
108: VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
109: (*ctx->f)(t,ctx->U,ctx->V,ctx->UF);
110: VecScatterBegin(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
111: VecScatterEnd(ctx->scatterU,ctx->UF,F,INSERT_VALUES,SCATTER_FORWARD);
112: return 0;
113: }
115: /*
116: Defines the nonlinear function that is passed to the time-integrator
118: */
119: PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
120: {
121: AppCtx *ctx = (AppCtx*)actx;
124: VecCopy(UVdot,F);
125: VecScatterBegin(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
126: VecScatterEnd(ctx->scatterU,UV,ctx->U,INSERT_VALUES,SCATTER_REVERSE);
127: VecScatterBegin(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
128: VecScatterEnd(ctx->scatterV,UV,ctx->V,INSERT_VALUES,SCATTER_REVERSE);
129: (*ctx->F)(t,ctx->U,ctx->V,ctx->VF);
130: VecScatterBegin(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
131: VecScatterEnd(ctx->scatterV,ctx->VF,F,INSERT_VALUES,SCATTER_FORWARD);
132: return 0;
133: }