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: {
 24:   PetscFunctionBeginUser;
 25:   PetscCall(PetscNew(tsdae));
 26:   (*tsdae)->comm = comm;
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

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

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

 50: PetscErrorCode TSDAESimpleDestroy(TSDAESimple *tsdae)
 51: {
 52:   PetscFunctionBeginUser;
 53:   PetscCall((*(*tsdae)->destroy)(*tsdae));
 54:   PetscCall(VecDestroy(&(*tsdae)->U));
 55:   PetscCall(VecDestroy(&(*tsdae)->V));
 56:   PetscCall(PetscFree(*tsdae));
 57:   PetscFunctionReturn(PETSC_SUCCESS);
 58: }

 60: PetscErrorCode TSDAESimpleSolve(TSDAESimple tsdae, Vec Usolution)
 61: {
 62:   PetscFunctionBeginUser;
 63:   PetscCall((*tsdae->solve)(tsdae, Usolution));
 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: PetscErrorCode TSDAESimpleSetFromOptions(TSDAESimple tsdae)
 68: {
 69:   PetscFunctionBeginUser;
 70:   PetscCall((*tsdae->setfromoptions)(PetscOptionsObject, tsdae));
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 98:   PetscFunctionBeginUser;
 99:   red->t = t;
100:   red->U = U;
101:   PetscCall(SNESSolve(red->snes, NULL, tsdae->V));
102:   PetscCall((*tsdae->f)(t, U, tsdae->V, F, tsdae->fctx));
103:   PetscFunctionReturn(PETSC_SUCCESS);
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;

115:   PetscFunctionBeginUser;
116:   PetscCall((*tsdae->F)(red->t, red->U, V, F, tsdae->Fctx));
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

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

124:   PetscFunctionBeginUser;
125:   PetscCall(TSSolve(red->ts, U));
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

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

133:   PetscFunctionBeginUser;
134:   PetscCall(TSSetFromOptions(red->ts));
135:   PetscCall(SNESSetFromOptions(red->snes));
136:   PetscFunctionReturn(PETSC_SUCCESS);
137: }

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

143:   PetscFunctionBeginUser;
144:   PetscCall(TSDestroy(&red->ts));
145:   PetscCall(SNESDestroy(&red->snes));
146:   PetscCall(PetscFree(red));
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

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

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

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

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

171:   PetscCall(SNESCreate(tsdae->comm, &red->snes));
172:   PetscCall(SNESSetOptionsPrefix(red->snes, "tsdaesimple_"));
173:   PetscCall(SNESSetFunction(red->snes, NULL, TSDAESimple_Reduced_SNESFunction, tsdae));
174:   PetscCall(SNESSetJacobian(red->snes, NULL, NULL, SNESComputeJacobianDefault, tsdae));
175:   PetscFunctionReturn(PETSC_SUCCESS);
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;

202:   PetscFunctionBeginUser;
203:   PetscCall(VecSet(F, 0.0));
204:   PetscCall(VecScatterBegin(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE));
205:   PetscCall(VecScatterEnd(full->scatterU, UV, tsdae->U, INSERT_VALUES, SCATTER_REVERSE));
206:   PetscCall(VecScatterBegin(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE));
207:   PetscCall(VecScatterEnd(full->scatterV, UV, tsdae->V, INSERT_VALUES, SCATTER_REVERSE));
208:   PetscCall((*tsdae->f)(t, tsdae->U, tsdae->V, full->UF, tsdae->fctx));
209:   PetscCall(VecScatterBegin(full->scatterU, full->UF, F, INSERT_VALUES, SCATTER_FORWARD));
210:   PetscCall(VecScatterEnd(full->scatterU, full->UF, F, INSERT_VALUES, SCATTER_FORWARD));
211:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

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

256:   PetscFunctionBeginUser;
257:   PetscCall(TSSetFromOptions(full->ts));
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

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

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

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

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

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

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

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

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

311:   PetscCall(TSSetRHSFunction(full->ts, tsrhs, TSDAESimple_Full_TSRHSFunction, tsdae));
312:   PetscCall(TSSetIFunction(full->ts, NULL, TSDAESimple_Full_TSIFunction, tsdae));
313:   PetscCall(VecDestroy(&tsrhs));
314:   PetscFunctionReturn(PETSC_SUCCESS);
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: {
325:   PetscFunctionBeginUser;
326:   PetscCall(VecWAXPY(F, 1.0, U, V));
327:   PetscFunctionReturn(PETSC_SUCCESS);
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: {
336:   PetscFunctionBeginUser;
337:   PetscCall(VecWAXPY(F, -1.0, V, U));
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

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

346:   PetscFunctionBeginUser;
347:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
348:   PetscCall(TSDAESimpleCreate(PETSC_COMM_WORLD, &tsdae));

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

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

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

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

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