Actual source code: ex10.c
petsc-3.6.1 2015-08-06
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)(PetscOptions*,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: };
24: PetscErrorCode TSDAESimpleCreate(MPI_Comm comm,TSDAESimple *tsdae)
25: {
29: PetscNew(tsdae);
30: (*tsdae)->comm = comm;
31: return(0);
32: }
36: PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
37: {
41: tsdae->f = f;
42: tsdae->U = U;
43: PetscObjectReference((PetscObject)U);
44: tsdae->fctx = ctx;
45: return(0);
46: }
50: PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
51: {
55: tsdae->F = F;
56: tsdae->V = V;
57: PetscObjectReference((PetscObject)V);
58: tsdae->Fctx = ctx;
59: return(0);
60: }
64: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
65: {
69: (*(*tsdae)->destroy)(*tsdae);
70: VecDestroy(&(*tsdae)->U);
71: VecDestroy(&(*tsdae)->V);
72: PetscFree(*tsdae);
73: return(0);
74: }
78: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution)
79: {
83: (*tsdae->solve)(tsdae,Usolution);
84: return(0);
85: }
89: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
90: {
94: (*tsdae->setfromoptions)(PetscOptionsObject,tsdae);
95: return(0);
96: }
98: /* ----------------------------------------------------------------------------*/
99: /*
100: Integrates the system by integrating the reduced ODE system and solving the
101: algebraic constraints at each stage with a separate SNES solve.
102: */
104: typedef struct {
105: PetscReal t;
106: TS ts;
107: SNES snes;
108: Vec U;
109: } TSDAESimple_Reduced;
113: /*
114: Defines the RHS function that is passed to the time-integrator.
116: Solves F(U,V) for V and then computes f(U,V)
118: */
119: PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
120: {
121: TSDAESimple tsdae = (TSDAESimple)actx;
122: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
123: PetscErrorCode ierr;
126: red->t = t;
127: red->U = U;
128: SNESSolve(red->snes,NULL,tsdae->V);
129: (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);
130: return(0);
131: }
135: /*
136: Defines the nonlinear function that is passed to the nonlinear solver
138: */
139: PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes,Vec V,Vec F,void *actx)
140: {
141: TSDAESimple tsdae = (TSDAESimple)actx;
142: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
143: PetscErrorCode ierr;
146: (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);
147: return(0);
148: }
153: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U)
154: {
155: PetscErrorCode ierr;
156: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
159: TSSolve(red->ts,U);
160: return(0);
161: }
165: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(PetscOptions *PetscOptionsObject,TSDAESimple tsdae)
166: {
167: PetscErrorCode ierr;
168: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
171: TSSetFromOptions(red->ts);
172: SNESSetFromOptions(red->snes);
173: return(0);
174: }
178: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
179: {
180: PetscErrorCode ierr;
181: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
184: TSDestroy(&red->ts);
185: SNESDestroy(&red->snes);
186: PetscFree(red);
187: return(0);
188: }
192: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
193: {
194: PetscErrorCode ierr;
195: TSDAESimple_Reduced *red;
196: Vec tsrhs;
199: PetscNew(&red);
200: tsdae->data = red;
202: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
203: tsdae->solve = TSDAESimpleSolve_Reduced;
204: tsdae->destroy = TSDAESimpleDestroy_Reduced;
206: TSCreate(tsdae->comm,&red->ts);
207: TSSetProblemType(red->ts,TS_NONLINEAR);
208: TSSetType(red->ts,TSEULER);
209: VecDuplicate(tsdae->U,&tsrhs);
210: TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);
211: VecDestroy(&tsrhs);
213: SNESCreate(tsdae->comm,&red->snes);
214: SNESSetOptionsPrefix(red->snes,"tsdaesimple_");
215: SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);
216: SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);
217: return(0);
218: }
221: /* ----------------------------------------------------------------------------*/
223: /*
224: Integrates the system by integrating directly the entire DAE system
225: */
227: typedef struct {
228: TS ts;
229: Vec UV,UF,VF;
230: VecScatter scatterU,scatterV;
231: } TSDAESimple_Full;
235: /*
236: Defines the RHS function that is passed to the time-integrator.
238: f(U,V)
239: 0
241: */
242: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
243: {
244: TSDAESimple tsdae = (TSDAESimple)actx;
245: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
246: PetscErrorCode ierr;
249: VecSet(F,0.0);
250: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
251: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
252: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
253: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
254: (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);
255: VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
256: VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
257: return(0);
258: }
262: /*
263: Defines the nonlinear function that is passed to the nonlinear solver
265: \dot{U}
266: F(U,V)
268: */
269: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
270: {
271: TSDAESimple tsdae = (TSDAESimple)actx;
272: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
273: PetscErrorCode ierr;
276: VecCopy(UVdot,F);
277: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
278: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
279: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
280: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
281: (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);
282: VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
283: VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
284: return(0);
285: }
290: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
291: {
292: PetscErrorCode ierr;
293: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
296: VecSet(full->UV,1.0);
297: VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
298: VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
299: TSSolve(full->ts,full->UV);
300: VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
301: VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
302: return(0);
303: }
307: PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptions *PetscOptionsObject,TSDAESimple tsdae)
308: {
309: PetscErrorCode ierr;
310: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
313: TSSetFromOptions(full->ts);
314: return(0);
315: }
319: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
320: {
321: PetscErrorCode ierr;
322: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
325: TSDestroy(&full->ts);
326: VecDestroy(&full->UV);
327: VecDestroy(&full->UF);
328: VecDestroy(&full->VF);
329: VecScatterDestroy(&full->scatterU);
330: VecScatterDestroy(&full->scatterV);
331: PetscFree(full);
332: return(0);
333: }
337: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
338: {
339: PetscErrorCode ierr;
340: TSDAESimple_Full *full;
341: Vec tsrhs;
342: PetscInt nU,nV,UVstart;
343: IS is;
346: PetscNew(&full);
347: tsdae->data = full;
349: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
350: tsdae->solve = TSDAESimpleSolve_Full;
351: tsdae->destroy = TSDAESimpleDestroy_Full;
353: TSCreate(tsdae->comm,&full->ts);
354: TSSetProblemType(full->ts,TS_NONLINEAR);
355: TSSetType(full->ts,TSROSW);
356: VecDuplicate(tsdae->U,&full->UF);
357: VecDuplicate(tsdae->V,&full->VF);
359: VecGetLocalSize(tsdae->U,&nU);
360: VecGetLocalSize(tsdae->V,&nV);
361: VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);
362: VecDuplicate(tsrhs,&full->UV);
364: VecGetOwnershipRange(tsrhs,&UVstart,NULL);
365: ISCreateStride(tsdae->comm,nU,UVstart,1,&is);
366: VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);
367: ISDestroy(&is);
368: ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);
369: VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);
370: ISDestroy(&is);
372: TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);
373: TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);
374: VecDestroy(&tsrhs);
375: return(0);
376: }
379: /* ----------------------------------------------------------------------------*/
384: /*
385: Simple example: f(U,V) = U + V
387: */
388: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
389: {
393: VecWAXPY(F,1.0,U,V);
394: return(0);
395: }
399: /*
400: Simple example: F(U,V) = U - V
402: */
403: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
404: {
408: VecWAXPY(F,-1.0,V,U);
409: return(0);
410: }
414: int main(int argc,char **argv)
415: {
417: TSDAESimple tsdae;
418: Vec U,V,Usolution;
420: PetscInitialize(&argc,&argv,(char*)0,help);
421: TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);
423: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
424: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);
425: TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);
426: TSDAESimpleSetIFunction(tsdae,V,F,NULL);
428: VecDuplicate(U,&Usolution);
429: VecSet(Usolution,1.0);
431: /* TSDAESimpleSetUp_Full(tsdae); */
432: TSDAESimpleSetUp_Reduced(tsdae);
434: TSDAESimpleSetFromOptions(tsdae);
435: TSDAESimpleSolve(tsdae,Usolution);
436: TSDAESimpleDestroy(&tsdae);
438: VecDestroy(&U);
439: VecDestroy(&Usolution);
440: VecDestroy(&V);
441: PetscFinalize();
442: return 0;
443: }