Actual source code: tsdiscgrad.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  1: /*
  2:   Code for timestepping with discrete gradient integrators
  3: */
  4: #include <petsc/private/tsimpl.h>
  5: #include <petscdm.h>

  7: PetscBool DGCite = PETSC_FALSE;
  8: const char DGCitation[] = "@article{Gonzalez1996,\n"
  9:                           "  title   = {Time integration and discrete Hamiltonian systems},\n"
 10:                           "  author  = {Oscar Gonzalez},\n"
 11:                           "  journal = {Journal of Nonlinear Science},\n"
 12:                           "  volume  = {6},\n"
 13:                           "  pages   = {449--467},\n"
 14:                           "  doi     = {10.1007/978-1-4612-1246-1_10},\n"
 15:                           "  year    = {1996}\n}\n";

 17: typedef struct {
 18:   PetscReal stage_time;
 19:   Vec       X0, X, Xdot;
 20:   PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *);
 21:   PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *);
 22:   PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *);
 23: } TS_DiscGrad;

 25: static PetscErrorCode TSDiscGradGetX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
 26: {
 27:   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;

 31:   if (X0) {
 32:     if (dm && dm != ts->dm) {DMGetNamedGlobalVector(dm, "TSDiscGrad_X0", X0);}
 33:     else                    {*X0  = ts->vec_sol;}
 34:   }
 35:   if (Xdot) {
 36:     if (dm && dm != ts->dm) {DMGetNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);}
 37:     else                    {*Xdot = dg->Xdot;}
 38:   }
 39:   return(0);
 40: }

 42: static PetscErrorCode TSDiscGradRestoreX0AndXdot(TS ts, DM dm, Vec *X0, Vec *Xdot)
 43: {

 47:   if (X0) {
 48:     if (dm && dm != ts->dm) {DMRestoreNamedGlobalVector(dm, "TSDiscGrad_X0", X0);}
 49:   }
 50:   if (Xdot) {
 51:     if (dm && dm != ts->dm) {DMRestoreNamedGlobalVector(dm, "TSDiscGrad_Xdot", Xdot);}
 52:   }
 53:   return(0);
 54: }

 56: static PetscErrorCode DMCoarsenHook_TSDiscGrad(DM fine, DM coarse, void *ctx)
 57: {
 59:   return(0);
 60: }

 62: static PetscErrorCode DMRestrictHook_TSDiscGrad(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, void *ctx)
 63: {
 64:   TS             ts = (TS) ctx;
 65:   Vec            X0, Xdot, X0_c, Xdot_c;

 69:   TSDiscGradGetX0AndXdot(ts, fine, &X0, &Xdot);
 70:   TSDiscGradGetX0AndXdot(ts, coarse, &X0_c, &Xdot_c);
 71:   MatRestrict(restrct, X0, X0_c);
 72:   MatRestrict(restrct, Xdot, Xdot_c);
 73:   VecPointwiseMult(X0_c, rscale, X0_c);
 74:   VecPointwiseMult(Xdot_c, rscale, Xdot_c);
 75:   TSDiscGradRestoreX0AndXdot(ts, fine, &X0, &Xdot);
 76:   TSDiscGradRestoreX0AndXdot(ts, coarse, &X0_c, &Xdot_c);
 77:   return(0);
 78: }

 80: static PetscErrorCode DMSubDomainHook_TSDiscGrad(DM dm, DM subdm, void *ctx)
 81: {
 83:   return(0);
 84: }

 86: static PetscErrorCode DMSubDomainRestrictHook_TSDiscGrad(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, void *ctx)
 87: {
 88:   TS             ts = (TS) ctx;
 89:   Vec            X0, Xdot, X0_sub, Xdot_sub;

 93:   TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);
 94:   TSDiscGradGetX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);

 96:   VecScatterBegin(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);
 97:   VecScatterEnd(gscat, X0, X0_sub, INSERT_VALUES, SCATTER_FORWARD);

 99:   VecScatterBegin(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);
100:   VecScatterEnd(gscat, Xdot, Xdot_sub, INSERT_VALUES, SCATTER_FORWARD);

102:   TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);
103:   TSDiscGradRestoreX0AndXdot(ts, subdm, &X0_sub, &Xdot_sub);
104:   return(0);
105: }

107: static PetscErrorCode TSSetUp_DiscGrad(TS ts)
108: {
109:   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;
110:   DM             dm;

114:   if (!dg->X)    {VecDuplicate(ts->vec_sol, &dg->X);}
115:   if (!dg->X0)   {VecDuplicate(ts->vec_sol, &dg->X0);}
116:   if (!dg->Xdot) {VecDuplicate(ts->vec_sol, &dg->Xdot);}

118:   TSGetDM(ts, &dm);
119:   DMCoarsenHookAdd(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);
120:   DMSubDomainHookAdd(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);
121:   return(0);
122: }

124: static PetscErrorCode TSSetFromOptions_DiscGrad(PetscOptionItems *PetscOptionsObject, TS ts)
125: {

129:   PetscOptionsHead(PetscOptionsObject, "Discrete Gradients ODE solver options");
130:   {
131:     //PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSDiscGradSetTheta",th->Theta,&th->Theta,NULL);
132:   }
133:   PetscOptionsTail();
134:   return(0);
135: }

137: static PetscErrorCode TSView_DiscGrad(TS ts,PetscViewer viewer)
138: {
139:   PetscBool      iascii;

143:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
144:   if (iascii) {
145:     PetscViewerASCIIPrintf(viewer,"  Discrete Gradients\n");
146:   }
147:   return(0);
148: }

150: static PetscErrorCode TSReset_DiscGrad(TS ts)
151: {
152:   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;

156:   VecDestroy(&dg->X);
157:   VecDestroy(&dg->X0);
158:   VecDestroy(&dg->Xdot);
159:   return(0);
160: }

162: static PetscErrorCode TSDestroy_DiscGrad(TS ts)
163: {
164:   DM             dm;

168:   TSReset_DiscGrad(ts);
169:   TSGetDM(ts, &dm);
170:   if (dm) {
171:     DMCoarsenHookRemove(dm, DMCoarsenHook_TSDiscGrad, DMRestrictHook_TSDiscGrad, ts);
172:     DMSubDomainHookRemove(dm, DMSubDomainHook_TSDiscGrad, DMSubDomainRestrictHook_TSDiscGrad, ts);
173:   }
174:   PetscFree(ts->data);
175:   PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",NULL);
176:   PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",NULL);
177:   return(0);
178: }

180: static PetscErrorCode TSInterpolate_DiscGrad(TS ts, PetscReal t, Vec X)
181: {
182:   TS_DiscGrad   *dg = (TS_DiscGrad*)ts->data;
183:   PetscReal      dt = t - ts->ptime;

187:   VecCopy(ts->vec_sol, dg->X);
188:   VecWAXPY(X, dt, dg->Xdot, dg->X);
189:   return(0);
190: }

192: static PetscErrorCode TSDiscGrad_SNESSolve(TS ts, Vec b, Vec x)
193: {
194:   SNES           snes;
195:   PetscInt       nits, lits;

199:   TSGetSNES(ts, &snes);
200:   SNESSolve(snes, b, x);
201:   SNESGetIterationNumber(snes, &nits);
202:   SNESGetLinearSolveIterations(snes, &lits);
203:   ts->snes_its += nits;
204:   ts->ksp_its  += lits;
205:   return(0);
206: }

208: static PetscErrorCode TSStep_DiscGrad(TS ts)
209: {
210:   TS_DiscGrad   *dg = (TS_DiscGrad *) ts->data;
211:   TSAdapt        adapt;
212:   TSStepStatus   status          = TS_STEP_INCOMPLETE;
213:   PetscInt       rejections      = 0;
214:   PetscBool      stageok, accept = PETSC_TRUE;
215:   PetscReal      next_time_step  = ts->time_step;

219:   TSGetAdapt(ts, &adapt);
220:   if (!ts->steprollback) {VecCopy(ts->vec_sol, dg->X0);}

222:   while (!ts->reason && status != TS_STEP_COMPLETE) {
223:     PetscReal shift = 1/(0.5*ts->time_step);

225:     dg->stage_time = ts->ptime + 0.5*ts->time_step;

227:     VecCopy(dg->X0, dg->X);
228:     TSPreStage(ts, dg->stage_time);
229:     TSDiscGrad_SNESSolve(ts, NULL, dg->X);
230:     TSPostStage(ts, dg->stage_time, 0, &dg->X);
231:     TSAdaptCheckStage(adapt, ts, dg->stage_time, dg->X, &stageok);
232:     if (!stageok) goto reject_step;

234:     status = TS_STEP_PENDING;
235:     VecAXPBYPCZ(dg->Xdot, -shift, shift, 0, dg->X0, dg->X);
236:     VecAXPY(ts->vec_sol, ts->time_step, dg->Xdot);
237:     TSAdaptChoose(adapt, ts, ts->time_step, NULL, &next_time_step, &accept);
238:     status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
239:     if (!accept) {
240:       VecCopy(dg->X0, ts->vec_sol);
241:       ts->time_step = next_time_step;
242:       goto reject_step;
243:     }
244:     ts->ptime    += ts->time_step;
245:     ts->time_step = next_time_step;
246:     break;

248:   reject_step:
249:     ts->reject++; accept = PETSC_FALSE;
250:     if (!ts->reason && ts->max_reject >= 0 && ++rejections > ts->max_reject) {
251:       ts->reason = TS_DIVERGED_STEP_REJECTED;
252:       PetscInfo2(ts, "Step=%D, step rejections %D greater than current TS allowed, stopping solve\n", ts->steps, rejections);
253:     }
254:   }
255:   return(0);
256: }

258: static PetscErrorCode TSGetStages_DiscGrad(TS ts, PetscInt *ns, Vec **Y)
259: {
260:   TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;

263:   if (ns) *ns = 1;
264:   if (Y)  *Y  = &(dg->X);
265:   return(0);
266: }

268: /*
269:   This defines the nonlinear equation that is to be solved with SNES
270:     G(U) = F[t0 + 0.5*dt, U, (U-U0)/dt] = 0
271: */
272: static PetscErrorCode SNESTSFormFunction_DiscGrad(SNES snes, Vec x, Vec y, TS ts)
273: {
274:   TS_DiscGrad   *dg    = (TS_DiscGrad *) ts->data;
275:   PetscReal      shift = 1/(0.5*ts->time_step);
276:   Vec            X0, Xdot;
277:   DM             dm, dmsave;

281:   SNESGetDM(snes, &dm);
282:   TSDiscGradGetX0AndXdot(ts, dm, &X0, &Xdot);
283:   VecAXPBYPCZ(Xdot, -shift, shift, 0, X0, x);

285:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
286:   dmsave = ts->dm;
287:   ts->dm = dm;
288:   TSComputeIFunction(ts, dg->stage_time, x, Xdot, y, PETSC_FALSE);
289:   ts->dm = dmsave;
290:   TSDiscGradRestoreX0AndXdot(ts, dm, &X0, &Xdot);
291:   return(0);
292: }

294: static PetscErrorCode SNESTSFormJacobian_DiscGrad(SNES snes, Vec x, Mat A, Mat B, TS ts)
295: {
296:   TS_DiscGrad   *dg    = (TS_DiscGrad *) ts->data;
297:   PetscReal      shift = 1/(0.5*ts->time_step);
298:   Vec            Xdot;
299:   DM             dm,dmsave;

303:   SNESGetDM(snes, &dm);
304:   /* Xdot has already been computed in SNESTSFormFunction_DiscGrad (SNES guarantees this) */
305:   TSDiscGradGetX0AndXdot(ts, dm, NULL, &Xdot);

307:   dmsave = ts->dm;
308:   ts->dm = dm;
309:   TSComputeIJacobian(ts, dg->stage_time, x, Xdot, shift, A, B, PETSC_FALSE);
310:   ts->dm = dmsave;
311:   TSDiscGradRestoreX0AndXdot(ts, dm, NULL, &Xdot);
312:   return(0);
313: }

315: static PetscErrorCode TSDiscGradGetFormulation_DiscGrad(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *))
316: {
317:   TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;

320:   *Sfunc = dg->Sfunc;
321:   *Ffunc = dg->Ffunc;
322:   *Gfunc = dg->Gfunc;
323:   return(0);
324: }

326: static PetscErrorCode TSDiscGradSetFormulation_DiscGrad(TS ts, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *))
327: {
328:   TS_DiscGrad *dg = (TS_DiscGrad *) ts->data;

331:   dg->Sfunc = Sfunc;
332:   dg->Ffunc = Ffunc;
333:   dg->Gfunc = Gfunc;
334:   return(0);
335: }

337: /*MC
338:   TSDISCGRAD - ODE solver using the discrete gradients version of the implicit midpoint method

340:   Level: beginner

342:   Notes: This is the implicit midpoint rule, with an optional term that guarantees the discrete gradient property. This
343:   timestepper applies to systems of the form
344: $ u_t = S(u) grad F(u)
345:   where S(u) is a linear operator, and F is a functional of u.

347: .seealso: TSCreate(), TSSetType(), TS, TSTHETA, TSDiscGradSetFormulation()
348: M*/
349: PETSC_EXTERN PetscErrorCode TSCreate_DiscGrad(TS ts)
350: {
351:   TS_DiscGrad       *th;

355:   PetscCitationsRegister(DGCitation, &DGCite);
356:   ts->ops->reset           = TSReset_DiscGrad;
357:   ts->ops->destroy         = TSDestroy_DiscGrad;
358:   ts->ops->view            = TSView_DiscGrad;
359:   ts->ops->setfromoptions  = TSSetFromOptions_DiscGrad;
360:   ts->ops->setup           = TSSetUp_DiscGrad;
361:   ts->ops->step            = TSStep_DiscGrad;
362:   ts->ops->interpolate     = TSInterpolate_DiscGrad;
363:   ts->ops->getstages       = TSGetStages_DiscGrad;
364:   ts->ops->snesfunction    = SNESTSFormFunction_DiscGrad;
365:   ts->ops->snesjacobian    = SNESTSFormJacobian_DiscGrad;
366:   ts->default_adapt_type   = TSADAPTNONE;

368:   ts->usessnes = PETSC_TRUE;

370:   PetscNewLog(ts,&th);
371:   ts->data = (void*)th;
372:   PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradGetFormulation_C",TSDiscGradGetFormulation_DiscGrad);
373:   PetscObjectComposeFunction((PetscObject)ts,"TSDiscGradSetFormulation_C",TSDiscGradSetFormulation_DiscGrad);
374:   return(0);
375: }

377: /*@C
378:   TSDiscGradGetFormulation - Get the construction method for S, F, and grad F from the formulation u_t = S grad F

380:   Not Collective

382:   Input Parameter:
383: . ts - timestepping context

385:   Output Parameters:
386: + Sfunc - constructor for the S matrix from the formulation
387: . Ffunc - functional F from the formulation
388: - Gfunc - constructor for the gradient of F from the formulation


391:   Calling sequence of Sfunc:
392: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *)

394:   Calling sequence of Ffunc:
395: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *)

397:   Calling sequence of Gfunc:
398: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *)

400:   Level: intermediate

402: .seealso: TSDiscGradSetFormulation()
403: @*/
404: PetscErrorCode TSDiscGradGetFormulation(TS ts, PetscErrorCode (**Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (**Ffunc)(TS, PetscReal, Vec, PetscScalar *, void *), PetscErrorCode (**Gfunc)(TS, PetscReal, Vec, Vec, void *))
405: {

413:   PetscUseMethod(ts,"TSDiscGradGetFormulation_C",(TS,PetscErrorCode(**Sfunc)(TS,PetscReal,Vec,Mat,void*),PetscErrorCode(**Ffunc)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode(**Gfunc)(TS,PetscReal,Vec,Vec,void*)),(ts,Sfunc,Ffunc,Gfunc));
414:   return(0);
415: }

417: /*@C
418:   TSDiscGradSetFormulation - Set the construction method for S, F, and grad F from the formulation u_t = S(u) grad F(u)

420:   Not Collective

422:   Input Parameters:
423: + ts    - timestepping context
424: . Sfunc - constructor for the S matrix from the formulation
425: . Ffunc - functional F from the formulation
426: - Gfunc - constructor for the gradient of F from the formulation

428:   Calling sequence of Sfunc:
429: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Mat S, void *)

431:   Calling sequence of Ffunc:
432: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, PetscScalar *F, void *)

434:   Calling sequence of Gfunc:
435: $ PetscErrorCode func(TS ts, PetscReal time, Vec u, Vec G, void *)

437:   Level: Intermediate

439: .seealso: TSDiscGradGetFormulation()
440: @*/
441: PetscErrorCode TSDiscGradSetFormulation(TS ts, PetscErrorCode (*Sfunc)(TS, PetscReal, Vec, Mat, void *), PetscErrorCode (*Ffunc)(TS, PetscReal, Vec , PetscScalar *, void *), PetscErrorCode (*Gfunc)(TS, PetscReal, Vec, Vec, void *), void *ctx)
442: {

450:   PetscTryMethod(ts,"TSDiscGradSetFormulation_C",(TS,PetscErrorCode(*Sfunc)(TS,PetscReal,Vec,Mat,void*),PetscErrorCode(*Ffunc)(TS,PetscReal,Vec,PetscScalar*,void*),PetscErrorCode(*Gfunc)(TS,PetscReal,Vec,Vec,void*)),(ts,Sfunc,Ffunc,Gfunc));
451:   return(0);
452: }