Actual source code: ex7.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 except the user provided functions take input values as a single vector instead of two vectors
10: */
12: /*
13: f(U,V) = U + V
15: */
16: PetscErrorCode f(PetscReal t,Vec UV,Vec F)
17: {
18: const PetscScalar *u,*v;
19: PetscScalar *f;
20: PetscInt n,i;
23: VecGetLocalSize(UV,&n);
24: n = n/2;
25: VecGetArrayRead(UV,&u);
26: v = u + n;
27: VecGetArrayWrite(F,&f);
28: for (i=0; i<n; i++) f[i] = u[i] + v[i];
29: VecRestoreArrayRead(UV,&u);
30: VecRestoreArrayWrite(F,&f);
31: return 0;
32: }
34: /*
35: F(U,V) = U - V
36: */
37: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
38: {
39: const PetscScalar *u,*v;
40: PetscScalar *f;
41: PetscInt n,i;
44: VecGetLocalSize(UV,&n);
45: n = n/2;
46: VecGetArrayRead(UV,&u);
47: v = u + n;
48: VecGetArrayWrite(F,&f);
49: for (i=0; i<n; i++) f[i] = u[i] - v[i];
50: VecRestoreArrayRead(UV,&u);
51: VecRestoreArrayWrite(F,&f);
52: return 0;
53: }
55: typedef struct {
56: PetscReal t;
57: SNES snes;
58: Vec UV,V;
59: VecScatter scatterU,scatterV;
60: PetscErrorCode (*f)(PetscReal,Vec,Vec);
61: PetscErrorCode (*F)(PetscReal,Vec,Vec);
62: } AppCtx;
64: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
65: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);
67: int main(int argc,char **argv)
68: {
69: AppCtx ctx;
70: TS ts;
71: Vec tsrhs,U;
72: IS is;
73: PetscInt i;
74: PetscMPIInt rank;
76: PetscInitialize(&argc,&argv,(char*)0,help);
77: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
78: TSCreate(PETSC_COMM_WORLD,&ts);
79: TSSetProblemType(ts,TS_NONLINEAR);
80: TSSetType(ts,TSEULER);
81: TSSetFromOptions(ts);
82: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&tsrhs);
83: VecDuplicate(tsrhs,&U);
84: TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
85: TSSetMaxTime(ts,1.0);
86: ctx.f = f;
88: SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
89: SNESSetFromOptions(ctx.snes);
90: SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
91: SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
92: ctx.F = F;
93: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);
95: /* Create scatters to move between separate U and V representation and UV representation of solution */
96: VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&ctx.UV);
97: i = 2*rank;
98: ISCreateGeneral(PETSC_COMM_WORLD,1,&i,PETSC_COPY_VALUES,&is);
99: VecScatterCreate(U,NULL,ctx.UV,is,&ctx.scatterU);
100: ISDestroy(&is);
101: i = 2*rank + 1;
102: ISCreateGeneral(PETSC_COMM_WORLD,1,&i,PETSC_COPY_VALUES,&is);
103: VecScatterCreate(ctx.V,NULL,ctx.UV,is,&ctx.scatterV);
104: ISDestroy(&is);
106: VecSet(U,1.0);
107: TSSolve(ts,U);
109: VecDestroy(&ctx.V);
110: VecDestroy(&ctx.UV);
111: VecScatterDestroy(&ctx.scatterU);
112: VecScatterDestroy(&ctx.scatterV);
113: VecDestroy(&tsrhs);
114: VecDestroy(&U);
115: SNESDestroy(&ctx.snes);
116: TSDestroy(&ts);
117: PetscFinalize();
118: return 0;
119: }
121: /*
122: Defines the RHS function that is passed to the time-integrator.
124: Solves F(U,V) for V and then computes f(U,V)
125: */
126: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
127: {
128: AppCtx *ctx = (AppCtx*)actx;
131: ctx->t = t;
132: VecScatterBegin(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
133: VecScatterEnd(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
134: SNESSolve(ctx->snes,NULL,ctx->V);
135: VecScatterBegin(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
136: VecScatterEnd(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
137: (*ctx->f)(t,ctx->UV,F);
138: return 0;
139: }
141: /*
142: Defines the nonlinear function that is passed to the nonlinear solver
144: */
145: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
146: {
147: AppCtx *ctx = (AppCtx*)actx;
150: VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
151: VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
152: (*ctx->F)(ctx->t,ctx->UV,F);
153: return 0;
154: }
156: /*TEST
158: test:
159: args: -ts_monitor -ts_view
161: TEST*/