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: PetscErrorCode ierr;
19: const PetscScalar *u,*v;
20: PetscScalar *f;
21: PetscInt n,i;
24: VecGetLocalSize(UV,&n);
25: n = n/2;
26: VecGetArrayRead(UV,&u);
27: v = u + n;
28: VecGetArrayWrite(F,&f);
29: for (i=0; i<n; i++) f[i] = u[i] + v[i];
30: VecRestoreArrayRead(UV,&u);
31: VecRestoreArrayWrite(F,&f);
32: return(0);
33: }
35: /*
36: F(U,V) = U - V
37: */
38: PetscErrorCode F(PetscReal t,Vec UV,Vec F)
39: {
40: PetscErrorCode ierr;
41: const PetscScalar *u,*v;
42: PetscScalar *f;
43: PetscInt n,i;
46: VecGetLocalSize(UV,&n);
47: n = n/2;
48: VecGetArrayRead(UV,&u);
49: v = u + n;
50: VecGetArrayWrite(F,&f);
51: for (i=0; i<n; i++) f[i] = u[i] - v[i];
52: VecRestoreArrayRead(UV,&u);
53: VecRestoreArrayWrite(F,&f);
54: return(0);
55: }
57: typedef struct {
58: PetscReal t;
59: SNES snes;
60: Vec UV,V;
61: VecScatter scatterU,scatterV;
62: PetscErrorCode (*f)(PetscReal,Vec,Vec);
63: PetscErrorCode (*F)(PetscReal,Vec,Vec);
64: } AppCtx;
66: extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*);
67: extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*);
69: int main(int argc,char **argv)
70: {
72: AppCtx ctx;
73: TS ts;
74: Vec tsrhs,U;
75: IS is;
76: PetscInt i;
77: PetscMPIInt rank;
79: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
80: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
81: TSCreate(PETSC_COMM_WORLD,&ts);
82: TSSetProblemType(ts,TS_NONLINEAR);
83: TSSetType(ts,TSEULER);
84: TSSetFromOptions(ts);
85: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&tsrhs);
86: VecDuplicate(tsrhs,&U);
87: TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx);
88: TSSetMaxTime(ts,1.0);
89: ctx.f = f;
91: SNESCreate(PETSC_COMM_WORLD,&ctx.snes);
92: SNESSetFromOptions(ctx.snes);
93: SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx);
94: SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx);
95: ctx.F = F;
96: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&ctx.V);
98: /* Create scatters to move between separate U and V representation and UV representation of solution */
99: VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&ctx.UV);
100: i = 2*rank;
101: ISCreateGeneral(PETSC_COMM_WORLD,1,&i,PETSC_COPY_VALUES,&is);
102: VecScatterCreate(U,NULL,ctx.UV,is,&ctx.scatterU);
103: ISDestroy(&is);
104: i = 2*rank + 1;
105: ISCreateGeneral(PETSC_COMM_WORLD,1,&i,PETSC_COPY_VALUES,&is);
106: VecScatterCreate(ctx.V,NULL,ctx.UV,is,&ctx.scatterV);
107: ISDestroy(&is);
109: VecSet(U,1.0);
110: TSSolve(ts,U);
112: VecDestroy(&ctx.V);
113: VecDestroy(&ctx.UV);
114: VecScatterDestroy(&ctx.scatterU);
115: VecScatterDestroy(&ctx.scatterV);
116: VecDestroy(&tsrhs);
117: VecDestroy(&U);
118: SNESDestroy(&ctx.snes);
119: TSDestroy(&ts);
120: PetscFinalize();
121: return ierr;
122: }
124: /*
125: Defines the RHS function that is passed to the time-integrator.
127: Solves F(U,V) for V and then computes f(U,V)
128: */
129: PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
130: {
131: AppCtx *ctx = (AppCtx*)actx;
135: ctx->t = t;
136: VecScatterBegin(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
137: VecScatterEnd(ctx->scatterU,U,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
138: SNESSolve(ctx->snes,NULL,ctx->V);
139: VecScatterBegin(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
140: VecScatterEnd(ctx->scatterV,ctx->V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
141: (*ctx->f)(t,ctx->UV,F);
142: return(0);
143: }
145: /*
146: Defines the nonlinear function that is passed to the nonlinear solver
148: */
149: PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx)
150: {
151: AppCtx *ctx = (AppCtx*)actx;
155: VecScatterBegin(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
156: VecScatterEnd(ctx->scatterV,V,ctx->UV,INSERT_VALUES,SCATTER_FORWARD);
157: (*ctx->F)(ctx->t,ctx->UV,F);
158: return(0);
159: }
161: /*TEST
163: test:
164: args: -ts_monitor -ts_view
166: TEST*/