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: {
 24:   PetscNew(tsdae);
 25:   (*tsdae)->comm = comm;
 26:   return 0;
 27: }

 29: PetscErrorCode TSDAESimpleSetRHSFunction(TSDAESimple tsdae,Vec U,PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
 30: {
 31:   tsdae->f    = f;
 32:   tsdae->U    = U;
 33:   PetscObjectReference((PetscObject)U);
 34:   tsdae->fctx = ctx;
 35:   return 0;
 36: }

 38: PetscErrorCode TSDAESimpleSetIFunction(TSDAESimple tsdae,Vec V,PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec,void*),void *ctx)
 39: {
 40:   tsdae->F    = F;
 41:   tsdae->V    = V;
 42:   PetscObjectReference((PetscObject)V);
 43:   tsdae->Fctx = ctx;
 44:   return 0;
 45: }

 47: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
 48: {
 49:   (*(*tsdae)->destroy)(*tsdae);
 50:   VecDestroy(&(*tsdae)->U);
 51:   VecDestroy(&(*tsdae)->V);
 52:   PetscFree(*tsdae);
 53:   return 0;
 54: }

 56: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae,Vec Usolution)
 57: {
 58:   (*tsdae->solve)(tsdae,Usolution);
 59:   return 0;
 60: }

 62: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
 63: {
 64:   (*tsdae->setfromoptions)(PetscOptionsObject,tsdae);
 65:   return 0;
 66: }

 68: /* ----------------------------------------------------------------------------*/
 69: /*
 70:       Integrates the system by integrating the reduced ODE system and solving the
 71:    algebraic constraints at each stage with a separate SNES solve.
 72: */

 74: typedef struct {
 75:   PetscReal t;
 76:   TS        ts;
 77:   SNES      snes;
 78:   Vec       U;
 79: } TSDAESimple_Reduced;

 81: /*
 82:    Defines the RHS function that is passed to the time-integrator.

 84:    Solves F(U,V) for V and then computes f(U,V)

 86: */
 87: PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx)
 88: {
 89:   TSDAESimple         tsdae = (TSDAESimple)actx;
 90:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

 93:   red->t = t;
 94:   red->U = U;
 95:   SNESSolve(red->snes,NULL,tsdae->V);
 96:   (*tsdae->f)(t,U,tsdae->V,F,tsdae->fctx);
 97:   return 0;
 98: }

100: /*
101:    Defines the nonlinear function that is passed to the nonlinear solver

103: */
104: PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes,Vec V,Vec F,void *actx)
105: {
106:   TSDAESimple         tsdae = (TSDAESimple)actx;
107:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

110:   (*tsdae->F)(red->t,red->U,V,F,tsdae->Fctx);
111:   return 0;
112: }

114: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae,Vec U)
115: {
116:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

118:   TSSolve(red->ts,U);
119:   return 0;
120: }

122: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
123: {
124:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

126:   TSSetFromOptions(red->ts);
127:   SNESSetFromOptions(red->snes);
128:   return 0;
129: }

131: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
132: {
133:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced*)tsdae->data;

135:   TSDestroy(&red->ts);
136:   SNESDestroy(&red->snes);
137:   PetscFree(red);
138:   return 0;
139: }

141: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
142: {
143:   TSDAESimple_Reduced *red;
144:   Vec                 tsrhs;

146:   PetscNew(&red);
147:   tsdae->data = red;

149:   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
150:   tsdae->solve          = TSDAESimpleSolve_Reduced;
151:   tsdae->destroy        = TSDAESimpleDestroy_Reduced;

153:   TSCreate(tsdae->comm,&red->ts);
154:   TSSetProblemType(red->ts,TS_NONLINEAR);
155:   TSSetType(red->ts,TSEULER);
156:   TSSetExactFinalTime(red->ts,TS_EXACTFINALTIME_STEPOVER);
157:   VecDuplicate(tsdae->U,&tsrhs);
158:   TSSetRHSFunction(red->ts,tsrhs,TSDAESimple_Reduced_TSFunction,tsdae);
159:   VecDestroy(&tsrhs);

161:   SNESCreate(tsdae->comm,&red->snes);
162:   SNESSetOptionsPrefix(red->snes,"tsdaesimple_");
163:   SNESSetFunction(red->snes,NULL,TSDAESimple_Reduced_SNESFunction,tsdae);
164:   SNESSetJacobian(red->snes,NULL,NULL,SNESComputeJacobianDefault,tsdae);
165:   return 0;
166: }

168: /* ----------------------------------------------------------------------------*/

170: /*
171:       Integrates the system by integrating directly the entire DAE system
172: */

174: typedef struct {
175:   TS         ts;
176:   Vec        UV,UF,VF;
177:   VecScatter scatterU,scatterV;
178: } TSDAESimple_Full;

180: /*
181:    Defines the RHS function that is passed to the time-integrator.

183:    f(U,V)
184:    0

186: */
187: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts,PetscReal t,Vec UV,Vec F,void *actx)
188: {
189:   TSDAESimple      tsdae = (TSDAESimple)actx;
190:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

193:   VecSet(F,0.0);
194:   VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
195:   VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
196:   VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
197:   VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
198:   (*tsdae->f)(t,tsdae->U,tsdae->V,full->UF,tsdae->fctx);
199:   VecScatterBegin(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
200:   VecScatterEnd(full->scatterU,full->UF,F,INSERT_VALUES,SCATTER_FORWARD);
201:   return 0;
202: }

204: /*
205:    Defines the nonlinear function that is passed to the nonlinear solver

207:    \dot{U}
208:    F(U,V)

210: */
211: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx)
212: {
213:   TSDAESimple       tsdae = (TSDAESimple)actx;
214:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

217:   VecCopy(UVdot,F);
218:   VecScatterBegin(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
219:   VecScatterEnd(full->scatterU,UV,tsdae->U,INSERT_VALUES,SCATTER_REVERSE);
220:   VecScatterBegin(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
221:   VecScatterEnd(full->scatterV,UV,tsdae->V,INSERT_VALUES,SCATTER_REVERSE);
222:   (*tsdae->F)(t,tsdae->U,tsdae->V,full->VF,tsdae->Fctx);
223:   VecScatterBegin(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
224:   VecScatterEnd(full->scatterV,full->VF,F,INSERT_VALUES,SCATTER_FORWARD);
225:   return 0;
226: }

228: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae,Vec U)
229: {
230:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

232:   VecSet(full->UV,1.0);
233:   VecScatterBegin(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
234:   VecScatterEnd(full->scatterU,U,full->UV,INSERT_VALUES,SCATTER_FORWARD);
235:   TSSolve(full->ts,full->UV);
236:   VecScatterBegin(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
237:   VecScatterEnd(full->scatterU,full->UV,U,INSERT_VALUES,SCATTER_REVERSE);
238:   return 0;
239: }

241: PetscErrorCode TSDAESimpleSetFromOptions_Full(PetscOptionItems *PetscOptionsObject,TSDAESimple tsdae)
242: {
243:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

245:   TSSetFromOptions(full->ts);
246:   return 0;
247: }

249: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
250: {
251:   TSDAESimple_Full *full = (TSDAESimple_Full*)tsdae->data;

253:   TSDestroy(&full->ts);
254:   VecDestroy(&full->UV);
255:   VecDestroy(&full->UF);
256:   VecDestroy(&full->VF);
257:   VecScatterDestroy(&full->scatterU);
258:   VecScatterDestroy(&full->scatterV);
259:   PetscFree(full);
260:   return 0;
261: }

263: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
264: {
265:   TSDAESimple_Full *full;
266:   Vec              tsrhs;
267:   PetscInt         nU,nV,UVstart;
268:   IS               is;

270:   PetscNew(&full);
271:   tsdae->data = full;

273:   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
274:   tsdae->solve          = TSDAESimpleSolve_Full;
275:   tsdae->destroy        = TSDAESimpleDestroy_Full;

277:   TSCreate(tsdae->comm,&full->ts);
278:   TSSetProblemType(full->ts,TS_NONLINEAR);
279:   TSSetType(full->ts,TSROSW);
280:   TSSetExactFinalTime(full->ts,TS_EXACTFINALTIME_STEPOVER);
281:   VecDuplicate(tsdae->U,&full->UF);
282:   VecDuplicate(tsdae->V,&full->VF);

284:   VecGetLocalSize(tsdae->U,&nU);
285:   VecGetLocalSize(tsdae->V,&nV);
286:   VecCreateMPI(tsdae->comm,nU+nV,PETSC_DETERMINE,&tsrhs);
287:   VecDuplicate(tsrhs,&full->UV);

289:   VecGetOwnershipRange(tsrhs,&UVstart,NULL);
290:   ISCreateStride(tsdae->comm,nU,UVstart,1,&is);
291:   VecScatterCreate(tsdae->U,NULL,tsrhs,is,&full->scatterU);
292:   ISDestroy(&is);
293:   ISCreateStride(tsdae->comm,nV,UVstart+nU,1,&is);
294:   VecScatterCreate(tsdae->V,NULL,tsrhs,is,&full->scatterV);
295:   ISDestroy(&is);

297:   TSSetRHSFunction(full->ts,tsrhs,TSDAESimple_Full_TSRHSFunction,tsdae);
298:   TSSetIFunction(full->ts,NULL,TSDAESimple_Full_TSIFunction,tsdae);
299:   VecDestroy(&tsrhs);
300:   return 0;
301: }

303: /* ----------------------------------------------------------------------------*/

305: /*
306:    Simple example:   f(U,V) = U + V

308: */
309: PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
310: {
312:   VecWAXPY(F,1.0,U,V);
313:   return 0;
314: }

316: /*
317:    Simple example: F(U,V) = U - V

319: */
320: PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F,void *ctx)
321: {
323:   VecWAXPY(F,-1.0,V,U);
324:   return 0;
325: }

327: int main(int argc,char **argv)
328: {
329:   TSDAESimple    tsdae;
330:   Vec            U,V,Usolution;

332:   PetscInitialize(&argc,&argv,(char*)0,help);
333:   TSDAESimpleCreate(PETSC_COMM_WORLD,&tsdae);

335:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&U);
336:   VecCreateMPI(PETSC_COMM_WORLD,1,PETSC_DETERMINE,&V);
337:   TSDAESimpleSetRHSFunction(tsdae,U,f,NULL);
338:   TSDAESimpleSetIFunction(tsdae,V,F,NULL);

340:   VecDuplicate(U,&Usolution);
341:   VecSet(Usolution,1.0);

343:   /*  TSDAESimpleSetUp_Full(tsdae); */
344:   TSDAESimpleSetUp_Reduced(tsdae);

346:   TSDAESimpleSetFromOptions(tsdae);
347:   TSDAESimpleSolve(tsdae,Usolution);
348:   TSDAESimpleDestroy(&tsdae);

350:   VecDestroy(&U);
351:   VecDestroy(&Usolution);
352:   VecDestroy(&V);
353:   PetscFinalize();
354:   return 0;
355: }