Actual source code: ex10.c
petsc-3.12.5 2020-03-29
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: }
135: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U)
136: {
137: PetscErrorCode ierr;
138: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
141: TSSolve(red->ts,U);
142: return(0);
143: }
145: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
146: {
147: PetscErrorCode ierr;
148: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
151: TSSetFromOptions(red->ts);
152: SNESSetFromOptions(red->snes);
153: return(0);
154: }
156: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
157: {
158: PetscErrorCode ierr;
159: TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
162: TSDestroy(&red->ts);
163: SNESDestroy(&red->snes);
164: PetscFree(red);
165: return(0);
166: }
168: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
169: {
170: PetscErrorCode ierr;
171: TSDAESimple_Reduced *red;
172: Vec tsrhs;
175: PetscNew(&red);
176: tsdae->data = red;
178: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
179: tsdae->solve = TSDAESimpleSolve_Reduced;
180: tsdae->destroy = TSDAESimpleDestroy_Reduced;
182: TSCreate(tsdae->comm,&red->ts);
183: TSSetProblemType(red->ts,TS_NONLINEAR);
184: TSSetType(red->ts,TSEULER);
185: TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER);
186: VecDuplicate(tsdae->U,&tsrhs);
187: TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);
188: VecDestroy(&tsrhs);
190: SNESCreate(tsdae->comm,&red->snes);
191: SNESSetOptionsPrefix(red->snes,"tsdaesimple_");
192: SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);
193: SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);
194: return(0);
195: }
198: /* ----------------------------------------------------------------------------*/
200: /*
201: Integrates the system by integrating directly the entire DAE system
202: */
204: typedef struct {
205: TS ts;
206: Vec UV,UF,VF;
207: VecScatter scatterU,scatterV;
208: } TSDAESimple_Full;
210: /*
211: Defines the RHS function that is passed to the time-integrator.
213: f(U,V)
214: 0
216: */
217: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
218: {
219: TSDAESimple tsdae = (TSDAESimple)actx;
220: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
221: PetscErrorCode ierr;
224: VecSet(F,0.0);
225: VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
226: VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
227: VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
228: VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
229: (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);
230: VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
231: VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
232: return(0);
233: }
235: /*
236: Defines the nonlinear function that is passed to the nonlinear solver
238: \dot{U}
239: F(U,V)
241: */
242: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
243: {
244: TSDAESimple tsdae = (TSDAESimple)actx;
245: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
246: PetscErrorCode ierr;
249: VecCopy(UVdot,F);
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->VF,tsdae->Fctx);
255: VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
256: VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
257: return(0);
258: }
261: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
262: {
263: PetscErrorCode ierr;
264: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
267: VecSet(full->UV,1.0);
268: VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
269: VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
270: TSSolve(full->ts,full->UV);
271: VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
272: VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
273: return(0);
274: }
276: PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
277: {
278: PetscErrorCode ierr;
279: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
282: TSSetFromOptions(full->ts);
283: return(0);
284: }
286: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
287: {
288: PetscErrorCode ierr;
289: TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
292: TSDestroy(&full->ts);
293: VecDestroy(&full->UV);
294: VecDestroy(&full->UF);
295: VecDestroy(&full->VF);
296: VecScatterDestroy(&full->scatterU);
297: VecScatterDestroy(&full->scatterV);
298: PetscFree(full);
299: return(0);
300: }
302: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
303: {
304: PetscErrorCode ierr;
305: TSDAESimple_Full *full;
306: Vec tsrhs;
307: PetscInt nU,nV,UVstart;
308: IS is;
311: PetscNew(&full);
312: tsdae->data = full;
314: tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
315: tsdae->solve = TSDAESimpleSolve_Full;
316: tsdae->destroy = TSDAESimpleDestroy_Full;
318: TSCreate(tsdae->comm,&full->ts);
319: TSSetProblemType(full->ts,TS_NONLINEAR);
320: TSSetType(full->ts,TSROSW);
321: TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER);
322: VecDuplicate(tsdae->U,&full->UF);
323: VecDuplicate(tsdae->V,&full->VF);
325: VecGetLocalSize(tsdae->U,&nU);
326: VecGetLocalSize(tsdae->V,&nV);
327: VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);
328: VecDuplicate(tsrhs,&full->UV);
330: VecGetOwnershipRange(tsrhs,&UVstart,NULL);
331: ISCreateStride(tsdae->comm,nU,UVstart,1,&is);
332: VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);
333: ISDestroy(&is);
334: ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);
335: VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);
336: ISDestroy(&is);
338: TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);
339: TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);
340: VecDestroy(&tsrhs);
341: return(0);
342: }
345: /* ----------------------------------------------------------------------------*/
348: /*
349: Simple example: f(U,V) = U + V
351: */
352: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
353: {
357: VecWAXPY(F,1.0,U,V);
358: return(0);
359: }
361: /*
362: Simple example: F(U,V) = U - V
364: */
365: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
366: {
370: VecWAXPY(F,-1.0,V,U);
371: return(0);
372: }
374: int main(int argc,char **argv)
375: {
377: TSDAESimple tsdae;
378: Vec U,V,Usolution;
380: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
381: TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);
383: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
384: VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);
385: TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);
386: TSDAESimpleSetIFunction(tsdae,V,F,NULL);
388: VecDuplicate(U,&Usolution);
389: VecSet(Usolution,1.0);
391: /* TSDAESimpleSetUp_Full(tsdae); */
392: TSDAESimpleSetUp_Reduced(tsdae);
394: TSDAESimpleSetFromOptions(tsdae);
395: TSDAESimpleSolve(tsdae,Usolution);
396: TSDAESimpleDestroy(&tsdae);
398: VecDestroy(&U);
399: VecDestroy(&Usolution);
400: VecDestroy(&V);
401: PetscFinalize();
402: return ierr;
403: }