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)(PetscOptionItems*,TSDAESimple);
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: {
24: PetscNew(tsdae);
25: (*tsdae)->comm = comm;
26: return 0;
27: }
29: PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
30: {
31: tsdae->f = f;
32: tsdae->U = U;
33: PetscObjectReference((PetscObject)U);
34: tsdae->fctx = ctx;
35: return 0;
36: }
38: PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
39: {
40: tsdae->F = F;
41: tsdae->V = V;
42: PetscObjectReference((PetscObject)V);
43: tsdae->Fctx = ctx;
44: return 0;
45: }
47: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
48: {
49: (*(*tsdae)->destroy)(*tsdae);
50: VecDestroy(&(*tsdae)->U);
51: VecDestroy(&(*tsdae)->V);
52: PetscFree(*tsdae);
53: return 0;
54: }
56: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution)
57: {
58: (*tsdae->solve)(tsdae,Usolution);
59: return 0;
60: }
62: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
63: {
64: (*tsdae->setfromoptions)(PetscOptionsObject,tsdae);
65: return 0;
66: }
68: /* ----------------------------------------------------------------------------*/
69: /*
70: Integrates the system by integrating the reduced ODE system and solving the
71: algebraic constraints at each stage with a separate SNES solve.
72: */
74: typedef struct {
75: PetscReal t;
76: TS ts;
77: SNES snes;
78: Vec U;
79: } TSDAESimple_Reduced;
81: /*
82: Defines the RHS function that is passed to the time-integrator.
84: Solves F(U,V) for V and then computes f(U,V)
86: */
87: PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
88: {
89: TSDAESimple tsdae = (TSDAESimple)actx;
90: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
93: red->t = t;
94: red->U = U;
95: SNESSolve(red->snes,NULL,tsdae->V);
96: (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);
97: return 0;
98: }
100: /*
101: Defines the nonlinear function that is passed to the nonlinear solver
103: */
104: PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes,Vec V,Vec F,void *actx)
105: {
106: TSDAESimple tsdae = (TSDAESimple)actx;
107: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
110: (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);
111: return 0;
112: }
114: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U)
115: {
116: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
118: TSSolve(red->ts,U);
119: return 0;
120: }
122: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
123: {
124: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
126: TSSetFromOptions(red->ts);
127: SNESSetFromOptions(red->snes);
128: return 0;
129: }
131: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
132: {
133: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
135: TSDestroy(&red->ts);
136: SNESDestroy(&red->snes);
137: PetscFree(red);
138: return 0;
139: }
141: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
142: {
143: TSDAESimple_Reduced *red;
144: Vec tsrhs;
146: PetscNew(&red);
147: tsdae->data = red;
149: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
150: tsdae->solve = TSDAESimpleSolve_Reduced;
151: tsdae->destroy = TSDAESimpleDestroy_Reduced;
153: TSCreate(tsdae->comm,&red->ts);
154: TSSetProblemType(red->ts,TS_NONLINEAR);
155: TSSetType(red->ts,TSEULER);
156: TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER);
157: VecDuplicate(tsdae->U,&tsrhs);
158: TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);
159: VecDestroy(&tsrhs);
161: SNESCreate(tsdae->comm,&red->snes);
162: SNESSetOptionsPrefix(red->snes,"tsdaesimple_");
163: SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);
164: SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);
165: return 0;
166: }
168: /* ----------------------------------------------------------------------------*/
170: /*
171: Integrates the system by integrating directly the entire DAE system
172: */
174: typedef struct {
175: TS ts;
176: Vec UV,UF,VF;
177: VecScatter scatterU,scatterV;
178: } TSDAESimple_Full;
180: /*
181: Defines the RHS function that is passed to the time-integrator.
183: f(U,V)
184: 0
186: */
187: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
188: {
189: TSDAESimple tsdae = (TSDAESimple)actx;
190: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
193: VecSet(F,0.0);
194: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
195: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
196: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
197: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
198: (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);
199: VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
200: VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
201: return 0;
202: }
204: /*
205: Defines the nonlinear function that is passed to the nonlinear solver
207: \dot{U}
208: F(U,V)
210: */
211: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
212: {
213: TSDAESimple tsdae = (TSDAESimple)actx;
214: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
217: VecCopy(UVdot,F);
218: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
219: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
220: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
221: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
222: (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);
223: VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
224: VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
225: return 0;
226: }
228: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
229: {
230: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
232: VecSet(full->UV,1.0);
233: VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
234: VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
235: TSSolve(full->ts,full->UV);
236: VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
237: VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
238: return 0;
239: }
241: PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
242: {
243: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
245: TSSetFromOptions(full->ts);
246: return 0;
247: }
249: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
250: {
251: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
253: TSDestroy(&full->ts);
254: VecDestroy(&full->UV);
255: VecDestroy(&full->UF);
256: VecDestroy(&full->VF);
257: VecScatterDestroy(&full->scatterU);
258: VecScatterDestroy(&full->scatterV);
259: PetscFree(full);
260: return 0;
261: }
263: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
264: {
265: TSDAESimple_Full *full;
266: Vec tsrhs;
267: PetscInt nU,nV,UVstart;
268: IS is;
270: PetscNew(&full);
271: tsdae->data = full;
273: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
274: tsdae->solve = TSDAESimpleSolve_Full;
275: tsdae->destroy = TSDAESimpleDestroy_Full;
277: TSCreate(tsdae->comm,&full->ts);
278: TSSetProblemType(full->ts,TS_NONLINEAR);
279: TSSetType(full->ts,TSROSW);
280: TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER);
281: VecDuplicate(tsdae->U,&full->UF);
282: VecDuplicate(tsdae->V,&full->VF);
284: VecGetLocalSize(tsdae->U,&nU);
285: VecGetLocalSize(tsdae->V,&nV);
286: VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);
287: VecDuplicate(tsrhs,&full->UV);
289: VecGetOwnershipRange(tsrhs,&UVstart,NULL);
290: ISCreateStride(tsdae->comm,nU,UVstart,1,&is);
291: VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);
292: ISDestroy(&is);
293: ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);
294: VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);
295: ISDestroy(&is);
297: TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);
298: TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);
299: VecDestroy(&tsrhs);
300: return 0;
301: }
303: /* ----------------------------------------------------------------------------*/
305: /*
306: Simple example: f(U,V) = U + V
308: */
309: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
310: {
312: VecWAXPY(F,1.0,U,V);
313: return 0;
314: }
316: /*
317: Simple example: F(U,V) = U - V
319: */
320: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
321: {
323: VecWAXPY(F,-1.0,V,U);
324: return 0;
325: }
327: int main(int argc,char **argv)
328: {
329: TSDAESimple tsdae;
330: Vec U,V,Usolution;
332: PetscInitialize(&argc,&argv,(char*)0,help);
333: TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);
335: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
336: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);
337: TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);
338: TSDAESimpleSetIFunction(tsdae,V,F,NULL);
340: VecDuplicate(U,&Usolution);
341: VecSet(Usolution,1.0);
343: /* TSDAESimpleSetUp_Full(tsdae); */
344: TSDAESimpleSetUp_Reduced(tsdae);
346: TSDAESimpleSetFromOptions(tsdae);
347: TSDAESimpleSolve(tsdae,Usolution);
348: TSDAESimpleDestroy(&tsdae);
350: VecDestroy(&U);
351: VecDestroy(&Usolution);
352: VecDestroy(&V);
353: PetscFinalize();
354: return 0;
355: }