xref: /petsc/src/ts/tests/ex10.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94) 
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";
4 
5 #include <petscts.h>
6 
7 /* ----------------------------------------------------------------------------*/
8 
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 };
21 
22 PetscErrorCode TSDAESimpleCreate(MPI_Comm comm,TSDAESimple *tsdae)
23 {
24   PetscErrorCode ierr;
25 
26   PetscFunctionBegin;
27   ierr           = PetscNew(tsdae);CHKERRQ(ierr);
28   (*tsdae)->comm = comm;
29   PetscFunctionReturn(0);
30 }
31 
32 PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
33 {
34   PetscErrorCode ierr;
35 
36   PetscFunctionBegin;
37   tsdae->f    = f;
38   tsdae->U    = U;
39   ierr        = PetscObjectReference((PetscObject)U);CHKERRQ(ierr);
40   tsdae->fctx = ctx;
41   PetscFunctionReturn(0);
42 }
43 
44 PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
45 {
46   PetscErrorCode ierr;
47 
48   PetscFunctionBegin;
49   tsdae->F    = F;
50   tsdae->V    = V;
51   ierr        = PetscObjectReference((PetscObject)V);CHKERRQ(ierr);
52   tsdae->Fctx = ctx;
53   PetscFunctionReturn(0);
54 }
55 
56 PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
57 {
58   PetscErrorCode ierr;
59 
60   PetscFunctionBegin;
61   ierr = (*(*tsdae)->destroy)(*tsdae);CHKERRQ(ierr);
62   ierr = VecDestroy(&(*tsdae)->U);CHKERRQ(ierr);
63   ierr = VecDestroy(&(*tsdae)->V);CHKERRQ(ierr);
64   ierr = PetscFree(*tsdae);CHKERRQ(ierr);
65   PetscFunctionReturn(0);
66 }
67 
68 PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution)
69 {
70   PetscErrorCode ierr;
71 
72   PetscFunctionBegin;
73   ierr = (*tsdae->solve)(tsdae,Usolution);CHKERRQ(ierr);
74   PetscFunctionReturn(0);
75 }
76 
77 PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
78 {
79   PetscErrorCode ierr;
80 
81   PetscFunctionBegin;
82   ierr = (*tsdae->setfromoptions)(PetscOptionsObject,tsdae);CHKERRQ(ierr);
83   PetscFunctionReturn(0);
84 }
85 
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 */
91 
92 typedef struct {
93   PetscReal t;
94   TS        ts;
95   SNES      snes;
96   Vec       U;
97 } TSDAESimple_Reduced;
98 
99 /*
100    Defines the RHS function that is passed to the time-integrator.
101 
102    Solves F(U,V) for V and then computes f(U,V)
103 
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;
110 
111   PetscFunctionBeginUser;
112   red->t = t;
113   red->U = U;
114   ierr   = SNESSolve(red->snes,NULL,tsdae->V);CHKERRQ(ierr);
115   ierr   = (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);CHKERRQ(ierr);
116   PetscFunctionReturn(0);
117 }
118 
119 /*
120    Defines the nonlinear function that is passed to the nonlinear solver
121 
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;
128 
129   PetscFunctionBeginUser;
130   ierr = (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);CHKERRQ(ierr);
131   PetscFunctionReturn(0);
132 }
133 
134 
135 PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U)
136 {
137   PetscErrorCode      ierr;
138   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
139 
140   PetscFunctionBegin;
141   ierr = TSSolve(red->ts,U);CHKERRQ(ierr);
142   PetscFunctionReturn(0);
143 }
144 
145 PetscErrorCode TSDAESimpleSetFromOptions_Reduced(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
146 {
147   PetscErrorCode      ierr;
148   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
149 
150   PetscFunctionBegin;
151   ierr = TSSetFromOptions(red->ts);CHKERRQ(ierr);
152   ierr = SNESSetFromOptions(red->snes);CHKERRQ(ierr);
153   PetscFunctionReturn(0);
154 }
155 
156 PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
157 {
158   PetscErrorCode      ierr;
159   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;
160 
161   PetscFunctionBegin;
162   ierr = TSDestroy(&red->ts);CHKERRQ(ierr);
163   ierr = SNESDestroy(&red->snes);CHKERRQ(ierr);
164   ierr = PetscFree(red);CHKERRQ(ierr);
165   PetscFunctionReturn(0);
166 }
167 
168 PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
169 {
170   PetscErrorCode      ierr;
171   TSDAESimple_Reduced *red;
172   Vec                 tsrhs;
173 
174   PetscFunctionBegin;
175   ierr = PetscNew(&red);CHKERRQ(ierr);
176   tsdae->data = red;
177 
178   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
179   tsdae->solve          = TSDAESimpleSolve_Reduced;
180   tsdae->destroy        = TSDAESimpleDestroy_Reduced;
181 
182   ierr = TSCreate(tsdae->comm,&red->ts);CHKERRQ(ierr);
183   ierr = TSSetProblemType(red->ts,TS_NONLINEAR);CHKERRQ(ierr);
184   ierr = TSSetType(red->ts,TSEULER);CHKERRQ(ierr);
185   ierr = TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
186   ierr = VecDuplicate(tsdae->U,&tsrhs);CHKERRQ(ierr);
187   ierr = TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);CHKERRQ(ierr);
188   ierr = VecDestroy(&tsrhs);CHKERRQ(ierr);
189 
190   ierr = SNESCreate(tsdae->comm,&red->snes);CHKERRQ(ierr);
191   ierr = SNESSetOptionsPrefix(red->snes,"tsdaesimple_");CHKERRQ(ierr);
192   ierr = SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);CHKERRQ(ierr);
193   ierr = SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);CHKERRQ(ierr);
194   PetscFunctionReturn(0);
195 }
196 
197 
198 /* ----------------------------------------------------------------------------*/
199 
200 /*
201       Integrates the system by integrating directly the entire DAE system
202 */
203 
204 typedef struct {
205   TS         ts;
206   Vec        UV,UF,VF;
207   VecScatter scatterU,scatterV;
208 } TSDAESimple_Full;
209 
210 /*
211    Defines the RHS function that is passed to the time-integrator.
212 
213    f(U,V)
214    0
215 
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;
222 
223   PetscFunctionBeginUser;
224   ierr = VecSet(F,0.0);CHKERRQ(ierr);
225   ierr = VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
226   ierr = VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
227   ierr = VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
228   ierr = VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
229   ierr = (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);CHKERRQ(ierr);
230   ierr = VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
231   ierr = VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
232   PetscFunctionReturn(0);
233 }
234 
235 /*
236    Defines the nonlinear function that is passed to the nonlinear solver
237 
238    \dot{U}
239    F(U,V)
240 
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;
247 
248   PetscFunctionBeginUser;
249   ierr = VecCopy(UVdot,F);CHKERRQ(ierr);
250   ierr = VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
251   ierr = VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
252   ierr = VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
253   ierr = VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
254   ierr = (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);CHKERRQ(ierr);
255   ierr = VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
256   ierr = VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
257   PetscFunctionReturn(0);
258 }
259 
260 
261 PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
262 {
263   PetscErrorCode   ierr;
264   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
265 
266   PetscFunctionBegin;
267   ierr = VecSet(full->UV,1.0);CHKERRQ(ierr);
268   ierr = VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
269   ierr = VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
270   ierr = TSSolve(full->ts,full->UV);CHKERRQ(ierr);
271   ierr = VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
272   ierr = VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
273   PetscFunctionReturn(0);
274 }
275 
276 PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
277 {
278   PetscErrorCode   ierr;
279   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
280 
281   PetscFunctionBegin;
282   ierr = TSSetFromOptions(full->ts);CHKERRQ(ierr);
283   PetscFunctionReturn(0);
284 }
285 
286 PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
287 {
288   PetscErrorCode   ierr;
289   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;
290 
291   PetscFunctionBegin;
292   ierr = TSDestroy(&full->ts);CHKERRQ(ierr);
293   ierr = VecDestroy(&full->UV);CHKERRQ(ierr);
294   ierr = VecDestroy(&full->UF);CHKERRQ(ierr);
295   ierr = VecDestroy(&full->VF);CHKERRQ(ierr);
296   ierr = VecScatterDestroy(&full->scatterU);CHKERRQ(ierr);
297   ierr = VecScatterDestroy(&full->scatterV);CHKERRQ(ierr);
298   ierr = PetscFree(full);CHKERRQ(ierr);
299   PetscFunctionReturn(0);
300 }
301 
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;
309 
310   PetscFunctionBegin;
311   ierr = PetscNew(&full);CHKERRQ(ierr);
312   tsdae->data = full;
313 
314   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
315   tsdae->solve          = TSDAESimpleSolve_Full;
316   tsdae->destroy        = TSDAESimpleDestroy_Full;
317 
318   ierr = TSCreate(tsdae->comm,&full->ts);CHKERRQ(ierr);
319   ierr = TSSetProblemType(full->ts,TS_NONLINEAR);CHKERRQ(ierr);
320   ierr = TSSetType(full->ts,TSROSW);CHKERRQ(ierr);
321   ierr = TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
322   ierr = VecDuplicate(tsdae->U,&full->UF);CHKERRQ(ierr);
323   ierr = VecDuplicate(tsdae->V,&full->VF);CHKERRQ(ierr);
324 
325   ierr = VecGetLocalSize(tsdae->U,&nU);CHKERRQ(ierr);
326   ierr = VecGetLocalSize(tsdae->V,&nV);CHKERRQ(ierr);
327   ierr = VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);CHKERRQ(ierr);
328   ierr = VecDuplicate(tsrhs,&full->UV);CHKERRQ(ierr);
329 
330   ierr = VecGetOwnershipRange(tsrhs,&UVstart,NULL);CHKERRQ(ierr);
331   ierr = ISCreateStride(tsdae->comm,nU,UVstart,1,&is);CHKERRQ(ierr);
332   ierr = VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);CHKERRQ(ierr);
333   ierr = ISDestroy(&is);CHKERRQ(ierr);
334   ierr = ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);CHKERRQ(ierr);
335   ierr = VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);CHKERRQ(ierr);
336   ierr = ISDestroy(&is);CHKERRQ(ierr);
337 
338   ierr = TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);CHKERRQ(ierr);
339   ierr = TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);CHKERRQ(ierr);
340   ierr = VecDestroy(&tsrhs);CHKERRQ(ierr);
341   PetscFunctionReturn(0);
342 }
343 
344 
345 /* ----------------------------------------------------------------------------*/
346 
347 
348 /*
349    Simple example:   f(U,V) = U + V
350 
351 */
352 PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
353 {
354   PetscErrorCode ierr;
355 
356   PetscFunctionBeginUser;
357   ierr = VecWAXPY(F,1.0,U,V);CHKERRQ(ierr);
358   PetscFunctionReturn(0);
359 }
360 
361 /*
362    Simple example: F(U,V) = U - V
363 
364 */
365 PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
366 {
367   PetscErrorCode ierr;
368 
369   PetscFunctionBeginUser;
370   ierr = VecWAXPY(F,-1.0,V,U);CHKERRQ(ierr);
371   PetscFunctionReturn(0);
372 }
373 
374 int main(int argc,char **argv)
375 {
376   PetscErrorCode ierr;
377   TSDAESimple    tsdae;
378   Vec            U,V,Usolution;
379 
380   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
381   ierr = TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);CHKERRQ(ierr);
382 
383   ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);CHKERRQ(ierr);
384   ierr = VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);CHKERRQ(ierr);
385   ierr = TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);CHKERRQ(ierr);
386   ierr = TSDAESimpleSetIFunction(tsdae,V,F,NULL);CHKERRQ(ierr);
387 
388   ierr = VecDuplicate(U,&Usolution);CHKERRQ(ierr);
389   ierr = VecSet(Usolution,1.0);CHKERRQ(ierr);
390 
391   /*  ierr = TSDAESimpleSetUp_Full(tsdae);CHKERRQ(ierr); */
392   ierr = TSDAESimpleSetUp_Reduced(tsdae);CHKERRQ(ierr);
393 
394   ierr = TSDAESimpleSetFromOptions(tsdae);CHKERRQ(ierr);
395   ierr = TSDAESimpleSolve(tsdae,Usolution);CHKERRQ(ierr);
396   ierr = TSDAESimpleDestroy(&tsdae);CHKERRQ(ierr);
397 
398   ierr = VecDestroy(&U);CHKERRQ(ierr);
399   ierr = VecDestroy(&Usolution);CHKERRQ(ierr);
400   ierr = VecDestroy(&V);CHKERRQ(ierr);
401   ierr = PetscFinalize();
402   return ierr;
403 }
404 
405 
406 
407