Actual source code: ex10.c
1: static char help[] = "Simple wrapper object to solve DAE of the form:\n\
2: \\dot{U} = f(U,V)\n\
3: F(U,V) = 0\n\n";
5: #include <petscts.h>
7: /* ----------------------------------------------------------------------------*/
9: typedef struct _p_TSDAESimple *TSDAESimple;
10: struct _p_TSDAESimple {
11: MPI_Comm comm;
12: PetscErrorCode (*setfromoptions)(TSDAESimple, PetscOptionItems *);
13: PetscErrorCode (*solve)(TSDAESimple, Vec);
14: PetscErrorCode (*destroy)(TSDAESimple);
15: Vec U, V;
16: PetscErrorCode (*f)(PetscReal, Vec, Vec, Vec, void *);
17: PetscErrorCode (*F)(PetscReal, Vec, Vec, Vec, void *);
18: void *fctx, *Fctx;
19: void *data;
20: };
22: PetscErrorCode TSDAESimpleCreate(MPI_Comm comm, TSDAESimple *tsdae)
23: {
25: PetscNew(tsdae);
26: (*tsdae)->comm = comm;
27: return 0;
28: }
30: PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae, Vec U, PetscErrorCode (*f)(PetscReal, Vec, Vec, Vec, void *), void *ctx)
31: {
33: tsdae->f = f;
34: tsdae->U = U;
35: PetscObjectReference((PetscObject)U);
36: tsdae->fctx = ctx;
37: return 0;
38: }
40: PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae, Vec V, PetscErrorCode (*F)(PetscReal, Vec, Vec, Vec, void *), void *ctx)
41: {
43: tsdae->F = F;
44: tsdae->V = V;
45: PetscObjectReference((PetscObject)V);
46: tsdae->Fctx = ctx;
47: return 0;
48: }
50: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
51: {
53: (*(*tsdae)->destroy)(*tsdae);
54: VecDestroy(&(*tsdae)->U);
55: VecDestroy(&(*tsdae)->V);
56: PetscFree(*tsdae);
57: return 0;
58: }
60: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae, Vec Usolution)
61: {
63: (*tsdae->solve)(tsdae, Usolution);
64: return 0;
65: }
67: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
68: {
70: (*tsdae->setfromoptions)(PetscOptionsObject, tsdae);
71: return 0;
72: }
74: /* ----------------------------------------------------------------------------*/
75: /*
76: Integrates the system by integrating the reduced ODE system and solving the
77: algebraic constraints at each stage with a separate SNES solve.
78: */
80: typedef struct {
81: PetscReal t;
82: TS ts;
83: SNES snes;
84: Vec U;
85: } TSDAESimple_Reduced;
87: /*
88: Defines the RHS function that is passed to the time-integrator.
90: Solves F(U,V) for V and then computes f(U,V)
92: */
93: PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts, PetscReal t, Vec U, Vec F, void *actx)
94: {
95: TSDAESimple tsdae = (TSDAESimple)actx;
96: TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;
99: red->t = t;
100: red->U = U;
101: SNESSolve(red->snes, NULL, tsdae->V);
102: (*tsdae->f)(t, U, tsdae->V, F, tsdae->fctx);
103: return 0;
104: }
106: /*
107: Defines the nonlinear function that is passed to the nonlinear solver
109: */
110: PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes, Vec V, Vec F, void *actx)
111: {
112: TSDAESimple tsdae = (TSDAESimple)actx;
113: TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;
116: (*tsdae->F)(red->t, red->U, V, F, tsdae->Fctx);
117: return 0;
118: }
120: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae, Vec U)
121: {
122: TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;
125: TSSolve(red->ts, U);
126: return 0;
127: }
129: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(TSDAESimple tsdae, PetscOptionItems *PetscOptionsObject)
130: {
131: TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;
134: TSSetFromOptions(red->ts);
135: SNESSetFromOptions(red->snes);
136: return 0;
137: }
139: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
140: {
141: TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;
144: TSDestroy(&red->ts);
145: SNESDestroy(&red->snes);
146: PetscFree(red);
147: return 0;
148: }
150: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
151: {
152: TSDAESimple_Reduced *red;
153: Vec tsrhs;
156: PetscNew(&red);
157: tsdae->data = red;
159: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
160: tsdae->solve = TSDAESimpleSolve_Reduced;
161: tsdae->destroy = TSDAESimpleDestroy_Reduced;
163: TSCreate(tsdae->comm, &red->ts);
164: TSSetProblemType(red->ts, TS_NONLINEAR);
165: TSSetType(red->ts, TSEULER);
166: TSSetExactFinalTime(red->ts, TS_EXACTFINALTIME_STEPOVER);
167: VecDuplicate(tsdae->U, &tsrhs);
168: TSSetRHSFunction(red->ts, tsrhs, TSDAESimple_Reduced_TSFunction, tsdae);
169: VecDestroy(&tsrhs);
171: SNESCreate(tsdae->comm, &red->snes);
172: SNESSetOptionsPrefix(red->snes, "tsdaesimple_");
173: SNESSetFunction(red->snes, NULL, TSDAESimple_Reduced_SNESFunction, tsdae);
174: SNESSetJacobian(red->snes, NULL, NULL, SNESComputeJacobianDefault, tsdae);
175: return 0;
176: }
178: /* ----------------------------------------------------------------------------*/
180: /*
181: Integrates the system by integrating directly the entire DAE system
182: */
184: typedef struct {
185: TS ts;
186: Vec UV, UF, VF;
187: VecScatter scatterU, scatterV;
188: } TSDAESimple_Full;
190: /*
191: Defines the RHS function that is passed to the time-integrator.
193: f(U,V)
194: 0
196: */
197: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts, PetscReal t, Vec UV, Vec F, void *actx)
198: {
199: TSDAESimple tsdae = (TSDAESimple)actx;
200: TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;
203: VecSet(F, 0.0);
204: VecScatterBegin(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
205: VecScatterEnd(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
206: VecScatterBegin(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
207: VecScatterEnd(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
208: (*tsdae->f)(t, tsdae->U, tsdae->V, full->UF, tsdae->fctx);
209: VecScatterBegin(full->scatterU, full->UF, F, INSERT_VALUES, SCATTER_FORWARD);
210: VecScatterEnd(full->scatterU, full->UF, F, INSERT_VALUES, SCATTER_FORWARD);
211: return 0;
212: }
214: /*
215: Defines the nonlinear function that is passed to the nonlinear solver
217: \dot{U}
218: F(U,V)
220: */
221: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx)
222: {
223: TSDAESimple tsdae = (TSDAESimple)actx;
224: TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;
227: VecCopy(UVdot, F);
228: VecScatterBegin(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
229: VecScatterEnd(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
230: VecScatterBegin(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
231: VecScatterEnd(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
232: (*tsdae->F)(t, tsdae->U, tsdae->V, full->VF, tsdae->Fctx);
233: VecScatterBegin(full->scatterV, full->VF, F, INSERT_VALUES, SCATTER_FORWARD);
234: VecScatterEnd(full->scatterV, full->VF, F, INSERT_VALUES, SCATTER_FORWARD);
235: return 0;
236: }
238: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae, Vec U)
239: {
240: TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;
243: VecSet(full->UV, 1.0);
244: VecScatterBegin(full->scatterU, U, full->UV, INSERT_VALUES, SCATTER_FORWARD);
245: VecScatterEnd(full->scatterU, U, full->UV, INSERT_VALUES, SCATTER_FORWARD);
246: TSSolve(full->ts, full->UV);
247: VecScatterBegin(full->scatterU, full->UV, U, INSERT_VALUES, SCATTER_REVERSE);
248: VecScatterEnd(full->scatterU, full->UV, U, INSERT_VALUES, SCATTER_REVERSE);
249: return 0;
250: }
252: PetscErrorCode TSDAESimpleSetFromOptions_Full(TSDAESimple tsdae, PetscOptionItems *PetscOptionsObject)
253: {
254: TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;
257: TSSetFromOptions(full->ts);
258: return 0;
259: }
261: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
262: {
263: TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;
266: TSDestroy(&full->ts);
267: VecDestroy(&full->UV);
268: VecDestroy(&full->UF);
269: VecDestroy(&full->VF);
270: VecScatterDestroy(&full->scatterU);
271: VecScatterDestroy(&full->scatterV);
272: PetscFree(full);
273: return 0;
274: }
276: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
277: {
278: TSDAESimple_Full *full;
279: Vec tsrhs;
280: PetscInt nU, nV, UVstart;
281: IS is;
284: PetscNew(&full);
285: tsdae->data = full;
287: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
288: tsdae->solve = TSDAESimpleSolve_Full;
289: tsdae->destroy = TSDAESimpleDestroy_Full;
291: TSCreate(tsdae->comm, &full->ts);
292: TSSetProblemType(full->ts, TS_NONLINEAR);
293: TSSetType(full->ts, TSROSW);
294: TSSetExactFinalTime(full->ts, TS_EXACTFINALTIME_STEPOVER);
295: VecDuplicate(tsdae->U, &full->UF);
296: VecDuplicate(tsdae->V, &full->VF);
298: VecGetLocalSize(tsdae->U, &nU);
299: VecGetLocalSize(tsdae->V, &nV);
300: VecCreateMPI(tsdae->comm, nU + nV, PETSC_DETERMINE, &tsrhs);
301: VecDuplicate(tsrhs, &full->UV);
303: VecGetOwnershipRange(tsrhs, &UVstart, NULL);
304: ISCreateStride(tsdae->comm, nU, UVstart, 1, &is);
305: VecScatterCreate(tsdae->U, NULL, tsrhs, is, &full->scatterU);
306: ISDestroy(&is);
307: ISCreateStride(tsdae->comm, nV, UVstart + nU, 1, &is);
308: VecScatterCreate(tsdae->V, NULL, tsrhs, is, &full->scatterV);
309: ISDestroy(&is);
311: TSSetRHSFunction(full->ts, tsrhs, TSDAESimple_Full_TSRHSFunction, tsdae);
312: TSSetIFunction(full->ts, NULL, TSDAESimple_Full_TSIFunction, tsdae);
313: VecDestroy(&tsrhs);
314: return 0;
315: }
317: /* ----------------------------------------------------------------------------*/
319: /*
320: Simple example: f(U,V) = U + V
322: */
323: PetscErrorCode f(PetscReal t, Vec U, Vec V, Vec F, void *ctx)
324: {
326: VecWAXPY(F, 1.0, U, V);
327: return 0;
328: }
330: /*
331: Simple example: F(U,V) = U - V
333: */
334: PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F, void *ctx)
335: {
337: VecWAXPY(F, -1.0, V, U);
338: return 0;
339: }
341: int main(int argc, char **argv)
342: {
343: TSDAESimple tsdae;
344: Vec U, V, Usolution;
347: PetscInitialize(&argc, &argv, (char *)0, help);
348: TSDAESimpleCreate(PETSC_COMM_WORLD, &tsdae);
350: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &U);
351: VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &V);
352: TSDAESimpleSetRHSFunction(tsdae, U, f, NULL);
353: TSDAESimpleSetIFunction(tsdae, V, F, NULL);
355: VecDuplicate(U, &Usolution);
356: VecSet(Usolution, 1.0);
358: /* TSDAESimpleSetUp_Full(tsdae); */
359: TSDAESimpleSetUp_Reduced(tsdae);
361: TSDAESimpleSetFromOptions(tsdae);
362: TSDAESimpleSolve(tsdae, Usolution);
363: TSDAESimpleDestroy(&tsdae);
365: VecDestroy(&U);
366: VecDestroy(&Usolution);
367: VecDestroy(&V);
368: PetscFinalize();
369: return 0;
370: }