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: }