Actual source code: ex10.c
petsc-3.7.3 2016-08-01
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: };
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(PetscOptionItems *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: TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER);
210: VecDuplicate(tsdae->U,&tsrhs);
211: TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);
212: VecDestroy(&tsrhs);
214: SNESCreate(tsdae->comm,&red->snes);
215: SNESSetOptionsPrefix(red->snes,"tsdaesimple_");
216: SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);
217: SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);
218: return(0);
219: }
222: /* ----------------------------------------------------------------------------*/
224: /*
225: Integrates the system by integrating directly the entire DAE system
226: */
228: typedef struct {
229: TS ts;
230: Vec UV,UF,VF;
231: VecScatter scatterU,scatterV;
232: } TSDAESimple_Full;
236: /*
237: Defines the RHS function that is passed to the time-integrator.
239: f(U,V)
240: 0
242: */
243: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
244: {
245: TSDAESimple tsdae = (TSDAESimple)actx;
246: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
247: PetscErrorCode ierr;
250: VecSet(F,0.0);
251: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
252: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
253: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
254: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
255: (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);
256: VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
257: VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
258: return(0);
259: }
263: /*
264: Defines the nonlinear function that is passed to the nonlinear solver
266: \dot{U}
267: F(U,V)
269: */
270: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
271: {
272: TSDAESimple tsdae = (TSDAESimple)actx;
273: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
274: PetscErrorCode ierr;
277: VecCopy(UVdot,F);
278: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
279: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
280: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
281: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
282: (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);
283: VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
284: VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
285: return(0);
286: }
291: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
292: {
293: PetscErrorCode ierr;
294: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
297: VecSet(full->UV,1.0);
298: VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
299: VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
300: TSSolve(full->ts,full->UV);
301: VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
302: VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
303: return(0);
304: }
308: PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
309: {
310: PetscErrorCode ierr;
311: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
314: TSSetFromOptions(full->ts);
315: return(0);
316: }
320: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
321: {
322: PetscErrorCode ierr;
323: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
326: TSDestroy(&full->ts);
327: VecDestroy(&full->UV);
328: VecDestroy(&full->UF);
329: VecDestroy(&full->VF);
330: VecScatterDestroy(&full->scatterU);
331: VecScatterDestroy(&full->scatterV);
332: PetscFree(full);
333: return(0);
334: }
338: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
339: {
340: PetscErrorCode ierr;
341: TSDAESimple_Full *full;
342: Vec tsrhs;
343: PetscInt nU,nV,UVstart;
344: IS is;
347: PetscNew(&full);
348: tsdae->data = full;
350: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
351: tsdae->solve = TSDAESimpleSolve_Full;
352: tsdae->destroy = TSDAESimpleDestroy_Full;
354: TSCreate(tsdae->comm,&full->ts);
355: TSSetProblemType(full->ts,TS_NONLINEAR);
356: TSSetType(full->ts,TSROSW);
357: TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER);
358: VecDuplicate(tsdae->U,&full->UF);
359: VecDuplicate(tsdae->V,&full->VF);
361: VecGetLocalSize(tsdae->U,&nU);
362: VecGetLocalSize(tsdae->V,&nV);
363: VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);
364: VecDuplicate(tsrhs,&full->UV);
366: VecGetOwnershipRange(tsrhs,&UVstart,NULL);
367: ISCreateStride(tsdae->comm,nU,UVstart,1,&is);
368: VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);
369: ISDestroy(&is);
370: ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);
371: VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);
372: ISDestroy(&is);
374: TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);
375: TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);
376: VecDestroy(&tsrhs);
377: return(0);
378: }
381: /* ----------------------------------------------------------------------------*/
386: /*
387: Simple example: f(U,V) = U + V
389: */
390: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
391: {
395: VecWAXPY(F,1.0,U,V);
396: return(0);
397: }
401: /*
402: Simple example: F(U,V) = U - V
404: */
405: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
406: {
410: VecWAXPY(F,-1.0,V,U);
411: return(0);
412: }
416: int main(int argc,char **argv)
417: {
419: TSDAESimple tsdae;
420: Vec U,V,Usolution;
422: PetscInitialize(&argc,&argv,(char*)0,help);
423: TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);
425: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
426: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);
427: TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);
428: TSDAESimpleSetIFunction(tsdae,V,F,NULL);
430: VecDuplicate(U,&Usolution);
431: VecSet(Usolution,1.0);
433: /* TSDAESimpleSetUp_Full(tsdae); */
434: TSDAESimpleSetUp_Reduced(tsdae);
436: TSDAESimpleSetFromOptions(tsdae);
437: TSDAESimpleSolve(tsdae,Usolution);
438: TSDAESimpleDestroy(&tsdae);
440: VecDestroy(&U);
441: VecDestroy(&Usolution);
442: VecDestroy(&V);
443: PetscFinalize();
444: return 0;
445: }