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;
77: PetscInitialize(&argc, &argv, (char *)0, help);
78: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
79: TSCreate(PETSC_COMM_WORLD, &ts);
80: TSSetProblemType(ts, TS_NONLINEAR);
81: TSSetType(ts, TSEULER);
82: TSSetFromOptions(ts);
83: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &tsrhs);
84: VecDuplicate(tsrhs, &U);
85: TSSetRHSFunction(ts, tsrhs, TSFunction, &ctx);
86: TSSetMaxTime(ts, 1.0);
87: ctx.f = f;
89: SNESCreate(PETSC_COMM_WORLD, &ctx.snes);
90: SNESSetFromOptions(ctx.snes);
91: SNESSetFunction(ctx.snes, NULL, SNESFunction, &ctx);
92: SNESSetJacobian(ctx.snes, NULL, NULL, SNESComputeJacobianDefault, &ctx);
93: ctx.F = F;
94: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.V);
96: /* Create scatters to move between separate U and V representation and UV representation of solution */
97: VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &ctx.UV);
98: i = 2 * rank;
99: ISCreateGeneral(PETSC_COMM_WORLD, 1, &i, PETSC_COPY_VALUES, &is);
100: VecScatterCreate(U, NULL, ctx.UV, is, &ctx.scatterU);
101: ISDestroy(&is);
102: i = 2 * rank + 1;
103: ISCreateGeneral(PETSC_COMM_WORLD, 1, &i, PETSC_COPY_VALUES, &is);
104: VecScatterCreate(ctx.V, NULL, ctx.UV, is, &ctx.scatterV);
105: ISDestroy(&is);
107: VecSet(U, 1.0);
108: TSSolve(ts, U);
110: VecDestroy(&ctx.V);
111: VecDestroy(&ctx.UV);
112: VecScatterDestroy(&ctx.scatterU);
113: VecScatterDestroy(&ctx.scatterV);
114: VecDestroy(&tsrhs);
115: VecDestroy(&U);
116: SNESDestroy(&ctx.snes);
117: TSDestroy(&ts);
118: PetscFinalize();
119: return 0;
120: }
122: /*
123: Defines the RHS function that is passed to the time-integrator.
125: Solves F(U,V) for V and then computes f(U,V)
126: */
127: PetscErrorCode TSFunction(TS ts, PetscReal t, Vec U, Vec F, void *actx)
128: {
129: AppCtx *ctx = (AppCtx *)actx;
132: ctx->t = t;
133: VecScatterBegin(ctx->scatterU, U, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
134: VecScatterEnd(ctx->scatterU, U, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
135: SNESSolve(ctx->snes, NULL, ctx->V);
136: VecScatterBegin(ctx->scatterV, ctx->V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
137: VecScatterEnd(ctx->scatterV, ctx->V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
138: (*ctx->f)(t, ctx->UV, F);
139: return 0;
140: }
142: /*
143: Defines the nonlinear function that is passed to the nonlinear solver
145: */
146: PetscErrorCode SNESFunction(SNES snes, Vec V, Vec F, void *actx)
147: {
148: AppCtx *ctx = (AppCtx *)actx;
151: VecScatterBegin(ctx->scatterV, V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
152: VecScatterEnd(ctx->scatterV, V, ctx->UV, INSERT_VALUES, SCATTER_FORWARD);
153: (*ctx->F)(ctx->t, ctx->UV, F);
154: return 0;
155: }
157: /*TEST
159: test:
160: args: -ts_monitor -ts_view
162: TEST*/