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

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

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

 50: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
 51: {
 53:   (*(*tsdae)->destroy)(*tsdae);
 54:   VecDestroy(&(*tsdae)->U);
 55:   VecDestroy(&(*tsdae)->V);
 56:   PetscFree(*tsdae);
 57:   return 0;
 58: }

 60: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae, Vec Usolution)
 61: {
 63:   (*tsdae->solve)(tsdae, Usolution);
 64:   return 0;
 65: }

 67: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
 68: {
 70:   (*tsdae->setfromoptions)(PetscOptionsObject, tsdae);
 71:   return 0;
 72: }

 74: /* ----------------------------------------------------------------------------*/
 75: /*
 76:       Integrates the system by integrating the reduced ODE system and solving the
 77:    algebraic constraints at each stage with a separate SNES solve.
 78: */

 80: typedef struct {
 81:   PetscReal t;
 82:   TS        ts;
 83:   SNES      snes;
 84:   Vec       U;
 85: } TSDAESimple_Reduced;

 87: /*
 88:    Defines the RHS function that is passed to the time-integrator.

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

 92: */
 93: PetscErrorCode TSDAESimple_Reduced_TSFunction(TS ts, PetscReal t, Vec U, Vec F, void *actx)
 94: {
 95:   TSDAESimple          tsdae = (TSDAESimple)actx;
 96:   TSDAESimple_Reduced *red   = (TSDAESimple_Reduced *)tsdae->data;

 99:   red->t = t;
100:   red->U = U;
101:   SNESSolve(red->snes, NULL, tsdae->V);
102:   (*tsdae->f)(t, U, tsdae->V, F, tsdae->fctx);
103:   return 0;
104: }

106: /*
107:    Defines the nonlinear function that is passed to the nonlinear solver

109: */
110: PetscErrorCode TSDAESimple_Reduced_SNESFunction(SNES snes, Vec V, Vec F, void *actx)
111: {
112:   TSDAESimple          tsdae = (TSDAESimple)actx;
113:   TSDAESimple_Reduced *red   = (TSDAESimple_Reduced *)tsdae->data;

116:   (*tsdae->F)(red->t, red->U, V, F, tsdae->Fctx);
117:   return 0;
118: }

120: PetscErrorCode TSDAESimpleSolve_Reduced(TSDAESimple tsdae, Vec U)
121: {
122:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;

125:   TSSolve(red->ts, U);
126:   return 0;
127: }

129: PetscErrorCode TSDAESimpleSetFromOptions_Reduced(TSDAESimple tsdae, PetscOptionItems *PetscOptionsObject)
130: {
131:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;

134:   TSSetFromOptions(red->ts);
135:   SNESSetFromOptions(red->snes);
136:   return 0;
137: }

139: PetscErrorCode TSDAESimpleDestroy_Reduced(TSDAESimple tsdae)
140: {
141:   TSDAESimple_Reduced *red = (TSDAESimple_Reduced *)tsdae->data;

144:   TSDestroy(&red->ts);
145:   SNESDestroy(&red->snes);
146:   PetscFree(red);
147:   return 0;
148: }

150: PetscErrorCode TSDAESimpleSetUp_Reduced(TSDAESimple tsdae)
151: {
152:   TSDAESimple_Reduced *red;
153:   Vec                  tsrhs;

156:   PetscNew(&red);
157:   tsdae->data = red;

159:   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Reduced;
160:   tsdae->solve          = TSDAESimpleSolve_Reduced;
161:   tsdae->destroy        = TSDAESimpleDestroy_Reduced;

163:   TSCreate(tsdae->comm, &red->ts);
164:   TSSetProblemType(red->ts, TS_NONLINEAR);
165:   TSSetType(red->ts, TSEULER);
166:   TSSetExactFinalTime(red->ts, TS_EXACTFINALTIME_STEPOVER);
167:   VecDuplicate(tsdae->U, &tsrhs);
168:   TSSetRHSFunction(red->ts, tsrhs, TSDAESimple_Reduced_TSFunction, tsdae);
169:   VecDestroy(&tsrhs);

171:   SNESCreate(tsdae->comm, &red->snes);
172:   SNESSetOptionsPrefix(red->snes, "tsdaesimple_");
173:   SNESSetFunction(red->snes, NULL, TSDAESimple_Reduced_SNESFunction, tsdae);
174:   SNESSetJacobian(red->snes, NULL, NULL, SNESComputeJacobianDefault, tsdae);
175:   return 0;
176: }

178: /* ----------------------------------------------------------------------------*/

180: /*
181:       Integrates the system by integrating directly the entire DAE system
182: */

184: typedef struct {
185:   TS         ts;
186:   Vec        UV, UF, VF;
187:   VecScatter scatterU, scatterV;
188: } TSDAESimple_Full;

190: /*
191:    Defines the RHS function that is passed to the time-integrator.

193:    f(U,V)
194:    0

196: */
197: PetscErrorCode TSDAESimple_Full_TSRHSFunction(TS ts, PetscReal t, Vec UV, Vec F, void *actx)
198: {
199:   TSDAESimple       tsdae = (TSDAESimple)actx;
200:   TSDAESimple_Full *full  = (TSDAESimple_Full *)tsdae->data;

203:   VecSet(F, 0.0);
204:   VecScatterBegin(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
205:   VecScatterEnd(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
206:   VecScatterBegin(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
207:   VecScatterEnd(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
208:   (*tsdae->f)(t, tsdae->U, tsdae->V, full->UF, tsdae->fctx);
209:   VecScatterBegin(full->scatterU, full->UF, F, INSERT_VALUES, SCATTER_FORWARD);
210:   VecScatterEnd(full->scatterU, full->UF, F, INSERT_VALUES, SCATTER_FORWARD);
211:   return 0;
212: }

214: /*
215:    Defines the nonlinear function that is passed to the nonlinear solver

217:    \dot{U}
218:    F(U,V)

220: */
221: PetscErrorCode TSDAESimple_Full_TSIFunction(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx)
222: {
223:   TSDAESimple       tsdae = (TSDAESimple)actx;
224:   TSDAESimple_Full *full  = (TSDAESimple_Full *)tsdae->data;

227:   VecCopy(UVdot, F);
228:   VecScatterBegin(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
229:   VecScatterEnd(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE);
230:   VecScatterBegin(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
231:   VecScatterEnd(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE);
232:   (*tsdae->F)(t, tsdae->U, tsdae->V, full->VF, tsdae->Fctx);
233:   VecScatterBegin(full->scatterV, full->VF, F, INSERT_VALUES, SCATTER_FORWARD);
234:   VecScatterEnd(full->scatterV, full->VF, F, INSERT_VALUES, SCATTER_FORWARD);
235:   return 0;
236: }

238: PetscErrorCode TSDAESimpleSolve_Full(TSDAESimple tsdae, Vec U)
239: {
240:   TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;

243:   VecSet(full->UV, 1.0);
244:   VecScatterBegin(full->scatterU, U, full->UV, INSERT_VALUES, SCATTER_FORWARD);
245:   VecScatterEnd(full->scatterU, U, full->UV, INSERT_VALUES, SCATTER_FORWARD);
246:   TSSolve(full->ts, full->UV);
247:   VecScatterBegin(full->scatterU, full->UV, U, INSERT_VALUES, SCATTER_REVERSE);
248:   VecScatterEnd(full->scatterU, full->UV, U, INSERT_VALUES, SCATTER_REVERSE);
249:   return 0;
250: }

252: PetscErrorCode TSDAESimpleSetFromOptions_Full(TSDAESimple tsdae, PetscOptionItems *PetscOptionsObject)
253: {
254:   TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;

257:   TSSetFromOptions(full->ts);
258:   return 0;
259: }

261: PetscErrorCode TSDAESimpleDestroy_Full(TSDAESimple tsdae)
262: {
263:   TSDAESimple_Full *full = (TSDAESimple_Full *)tsdae->data;

266:   TSDestroy(&full->ts);
267:   VecDestroy(&full->UV);
268:   VecDestroy(&full->UF);
269:   VecDestroy(&full->VF);
270:   VecScatterDestroy(&full->scatterU);
271:   VecScatterDestroy(&full->scatterV);
272:   PetscFree(full);
273:   return 0;
274: }

276: PetscErrorCode TSDAESimpleSetUp_Full(TSDAESimple tsdae)
277: {
278:   TSDAESimple_Full *full;
279:   Vec               tsrhs;
280:   PetscInt          nU, nV, UVstart;
281:   IS                is;

284:   PetscNew(&full);
285:   tsdae->data = full;

287:   tsdae->setfromoptions = TSDAESimpleSetFromOptions_Full;
288:   tsdae->solve          = TSDAESimpleSolve_Full;
289:   tsdae->destroy        = TSDAESimpleDestroy_Full;

291:   TSCreate(tsdae->comm, &full->ts);
292:   TSSetProblemType(full->ts, TS_NONLINEAR);
293:   TSSetType(full->ts, TSROSW);
294:   TSSetExactFinalTime(full->ts, TS_EXACTFINALTIME_STEPOVER);
295:   VecDuplicate(tsdae->U, &full->UF);
296:   VecDuplicate(tsdae->V, &full->VF);

298:   VecGetLocalSize(tsdae->U, &nU);
299:   VecGetLocalSize(tsdae->V, &nV);
300:   VecCreateMPI(tsdae->comm, nU + nV, PETSC_DETERMINE, &tsrhs);
301:   VecDuplicate(tsrhs, &full->UV);

303:   VecGetOwnershipRange(tsrhs, &UVstart, NULL);
304:   ISCreateStride(tsdae->comm, nU, UVstart, 1, &is);
305:   VecScatterCreate(tsdae->U, NULL, tsrhs, is, &full->scatterU);
306:   ISDestroy(&is);
307:   ISCreateStride(tsdae->comm, nV, UVstart + nU, 1, &is);
308:   VecScatterCreate(tsdae->V, NULL, tsrhs, is, &full->scatterV);
309:   ISDestroy(&is);

311:   TSSetRHSFunction(full->ts, tsrhs, TSDAESimple_Full_TSRHSFunction, tsdae);
312:   TSSetIFunction(full->ts, NULL, TSDAESimple_Full_TSIFunction, tsdae);
313:   VecDestroy(&tsrhs);
314:   return 0;
315: }

317: /* ----------------------------------------------------------------------------*/

319: /*
320:    Simple example:   f(U,V) = U + V

322: */
323: PetscErrorCode f(PetscReal t, Vec U, Vec V, Vec F, void *ctx)
324: {
326:   VecWAXPY(F, 1.0, U, V);
327:   return 0;
328: }

330: /*
331:    Simple example: F(U,V) = U - V

333: */
334: PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F, void *ctx)
335: {
337:   VecWAXPY(F, -1.0, V, U);
338:   return 0;
339: }

341: int main(int argc, char **argv)
342: {
343:   TSDAESimple tsdae;
344:   Vec         U, V, Usolution;

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

350:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &U);
351:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &V);
352:   TSDAESimpleSetRHSFunction(tsdae, U, f, NULL);
353:   TSDAESimpleSetIFunction(tsdae, V, F, NULL);

355:   VecDuplicate(U, &Usolution);
356:   VecSet(Usolution, 1.0);

358:   /*  TSDAESimpleSetUp_Full(tsdae); */
359:   TSDAESimpleSetUp_Reduced(tsdae);

361:   TSDAESimpleSetFromOptions(tsdae);
362:   TSDAESimpleSolve(tsdae, Usolution);
363:   TSDAESimpleDestroy(&tsdae);

365:   VecDestroy(&U);
366:   VecDestroy(&Usolution);
367:   VecDestroy(&V);
368:   PetscFinalize();
369:   return 0;
370: }