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: {
27: PetscNew(tsdae);
28: (*tsdae)->comm = comm;
29: return(0);
30: }
32: PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
33: {
37: tsdae->f = f;
38: tsdae->U = U;
39: PetscObjectReference((PetscObject)U);
40: tsdae->fctx = ctx;
41: return(0);
42: }
44: PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
45: {
49: tsdae->F = F;
50: tsdae->V = V;
51: PetscObjectReference((PetscObject)V);
52: tsdae->Fctx = ctx;
53: return(0);
54: }
56: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
57: {
61: (*(*tsdae)->destroy)(*tsdae);
62: VecDestroy(&(*tsdae)->U);
63: VecDestroy(&(*tsdae)->V);
64: PetscFree(*tsdae);
65: return(0);
66: }
68: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution)
69: {
73: (*tsdae->solve)(tsdae,Usolution);
74: return(0);
75: }
77: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
78: {
82: (*tsdae->setfromoptions)(PetscOptionsObject,tsdae);
83: return(0);
84: }
86: /* ----------------------------------------------------------------------------*/
87: /*
88: Integrates the system by integrating the reduced ODE system and solving the
89: algebraic constraints at each stage with a separate SNES solve.
90: */
92: typedef struct {
93: PetscReal t;
94: TS ts;
95: SNES snes;
96: Vec U;
97: } TSDAESimple_Reduced;
99: /*
100: Defines the RHS function that is passed to the time-integrator.
102: Solves F(U,V) for V and then computes f(U,V)
104: */
105: PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
106: {
107: TSDAESimple tsdae = (TSDAESimple)actx;
108: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
109: PetscErrorCode ierr;
112: red->t = t;
113: red->U = U;
114: SNESSolve(red->snes,NULL,tsdae->V);
115: (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);
116: return(0);
117: }
119: /*
120: Defines the nonlinear function that is passed to the nonlinear solver
122: */
123: PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes,Vec V,Vec F,void *actx)
124: {
125: TSDAESimple tsdae = (TSDAESimple)actx;
126: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
127: PetscErrorCode ierr;
130: (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);
131: return(0);
132: }
134: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U)
135: {
136: PetscErrorCode ierr;
137: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
140: TSSolve(red->ts,U);
141: return(0);
142: }
144: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
145: {
146: PetscErrorCode ierr;
147: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
150: TSSetFromOptions(red->ts);
151: SNESSetFromOptions(red->snes);
152: return(0);
153: }
155: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
156: {
157: PetscErrorCode ierr;
158: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
161: TSDestroy(&red->ts);
162: SNESDestroy(&red->snes);
163: PetscFree(red);
164: return(0);
165: }
167: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
168: {
169: PetscErrorCode ierr;
170: TSDAESimple_Reduced *red;
171: Vec tsrhs;
174: PetscNew(&red);
175: tsdae->data = red;
177: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
178: tsdae->solve = TSDAESimpleSolve_Reduced;
179: tsdae->destroy = TSDAESimpleDestroy_Reduced;
181: TSCreate(tsdae->comm,&red->ts);
182: TSSetProblemType(red->ts,TS_NONLINEAR);
183: TSSetType(red->ts,TSEULER);
184: TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER);
185: VecDuplicate(tsdae->U,&tsrhs);
186: TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);
187: VecDestroy(&tsrhs);
189: SNESCreate(tsdae->comm,&red->snes);
190: SNESSetOptionsPrefix(red->snes,"tsdaesimple_");
191: SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);
192: SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);
193: return(0);
194: }
196: /* ----------------------------------------------------------------------------*/
198: /*
199: Integrates the system by integrating directly the entire DAE system
200: */
202: typedef struct {
203: TS ts;
204: Vec UV,UF,VF;
205: VecScatter scatterU,scatterV;
206: } TSDAESimple_Full;
208: /*
209: Defines the RHS function that is passed to the time-integrator.
211: f(U,V)
212: 0
214: */
215: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
216: {
217: TSDAESimple tsdae = (TSDAESimple)actx;
218: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
219: PetscErrorCode ierr;
222: VecSet(F,0.0);
223: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
224: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
225: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
226: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
227: (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);
228: VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
229: VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
230: return(0);
231: }
233: /*
234: Defines the nonlinear function that is passed to the nonlinear solver
236: \dot{U}
237: F(U,V)
239: */
240: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
241: {
242: TSDAESimple tsdae = (TSDAESimple)actx;
243: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
244: PetscErrorCode ierr;
247: VecCopy(UVdot,F);
248: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
249: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
250: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
251: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
252: (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);
253: VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
254: VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
255: return(0);
256: }
258: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
259: {
260: PetscErrorCode ierr;
261: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
264: VecSet(full->UV,1.0);
265: VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
266: VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
267: TSSolve(full->ts,full->UV);
268: VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
269: VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
270: return(0);
271: }
273: PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
274: {
275: PetscErrorCode ierr;
276: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
279: TSSetFromOptions(full->ts);
280: return(0);
281: }
283: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
284: {
285: PetscErrorCode ierr;
286: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
289: TSDestroy(&full->ts);
290: VecDestroy(&full->UV);
291: VecDestroy(&full->UF);
292: VecDestroy(&full->VF);
293: VecScatterDestroy(&full->scatterU);
294: VecScatterDestroy(&full->scatterV);
295: PetscFree(full);
296: return(0);
297: }
299: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
300: {
301: PetscErrorCode ierr;
302: TSDAESimple_Full *full;
303: Vec tsrhs;
304: PetscInt nU,nV,UVstart;
305: IS is;
308: PetscNew(&full);
309: tsdae->data = full;
311: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
312: tsdae->solve = TSDAESimpleSolve_Full;
313: tsdae->destroy = TSDAESimpleDestroy_Full;
315: TSCreate(tsdae->comm,&full->ts);
316: TSSetProblemType(full->ts,TS_NONLINEAR);
317: TSSetType(full->ts,TSROSW);
318: TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER);
319: VecDuplicate(tsdae->U,&full->UF);
320: VecDuplicate(tsdae->V,&full->VF);
322: VecGetLocalSize(tsdae->U,&nU);
323: VecGetLocalSize(tsdae->V,&nV);
324: VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);
325: VecDuplicate(tsrhs,&full->UV);
327: VecGetOwnershipRange(tsrhs,&UVstart,NULL);
328: ISCreateStride(tsdae->comm,nU,UVstart,1,&is);
329: VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);
330: ISDestroy(&is);
331: ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);
332: VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);
333: ISDestroy(&is);
335: TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);
336: TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);
337: VecDestroy(&tsrhs);
338: return(0);
339: }
341: /* ----------------------------------------------------------------------------*/
343: /*
344: Simple example: f(U,V) = U + V
346: */
347: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
348: {
352: VecWAXPY(F,1.0,U,V);
353: return(0);
354: }
356: /*
357: Simple example: F(U,V) = U - V
359: */
360: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
361: {
365: VecWAXPY(F,-1.0,V,U);
366: return(0);
367: }
369: int main(int argc,char **argv)
370: {
372: TSDAESimple tsdae;
373: Vec U,V,Usolution;
375: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
376: TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);
378: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
379: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);
380: TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);
381: TSDAESimpleSetIFunction(tsdae,V,F,NULL);
383: VecDuplicate(U,&Usolution);
384: VecSet(Usolution,1.0);
386: /* TSDAESimpleSetUp_Full(tsdae); */
387: TSDAESimpleSetUp_Reduced(tsdae);
389: TSDAESimpleSetFromOptions(tsdae);
390: TSDAESimpleSolve(tsdae,Usolution);
391: TSDAESimpleDestroy(&tsdae);
393: VecDestroy(&U);
394: VecDestroy(&Usolution);
395: VecDestroy(&V);
396: PetscFinalize();
397: return ierr;
398: }