Actual source code: theta.c

petsc-master 2016-02-06
Report Typos and Errors
  1: /*
  2:   Code for timestepping with implicit Theta method
  3: */
  4: #include <petsc/private/tsimpl.h>                /*I   "petscts.h"   I*/
  5: #include <petscsnes.h>
  6: #include <petscdm.h>
  7: #include <petscmat.h>

  9: typedef struct {
 10:   Vec          X,Xdot;                   /* Storage for one stage */
 11:   Vec          X0;                       /* work vector to store X0 */
 12:   Vec          affine;                   /* Affine vector needed for residual at beginning of step */
 13:   Vec          *VecsDeltaLam;             /* Increment of the adjoint sensitivity w.r.t IC at stage*/
 14:   Vec          *VecsDeltaMu;              /* Increment of the adjoint sensitivity w.r.t P at stage*/
 15:   Vec          *VecsSensiTemp;            /* Vector to be timed with Jacobian transpose*/
 16:   PetscBool    extrapolate;
 17:   PetscBool    endpoint;
 18:   PetscReal    Theta;
 19:   PetscReal    stage_time;
 20:   TSStepStatus status;
 21:   char         *name;
 22:   PetscInt     order;
 23:   PetscReal    ccfl;               /* Placeholder for CFL coefficient relative to forward Euler */
 24:   PetscBool    adapt;  /* use time-step adaptivity ? */
 25:   PetscReal    ptime;
 26: } TS_Theta;

 30: static PetscErrorCode TSThetaGetX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
 31: {
 32:   TS_Theta       *th = (TS_Theta*)ts->data;

 36:   if (X0) {
 37:     if (dm && dm != ts->dm) {
 38:       DMGetNamedGlobalVector(dm,"TSTheta_X0",X0);
 39:     } else *X0 = ts->vec_sol;
 40:   }
 41:   if (Xdot) {
 42:     if (dm && dm != ts->dm) {
 43:       DMGetNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
 44:     } else *Xdot = th->Xdot;
 45:   }
 46:   return(0);
 47: }

 51: static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
 52: {

 56:   if (X0) {
 57:     if (dm && dm != ts->dm) {
 58:       DMRestoreNamedGlobalVector(dm,"TSTheta_X0",X0);
 59:     }
 60:   }
 61:   if (Xdot) {
 62:     if (dm && dm != ts->dm) {
 63:       DMRestoreNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
 64:     }
 65:   }
 66:   return(0);
 67: }

 71: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine,DM coarse,void *ctx)
 72: {

 75:   return(0);
 76: }

 80: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
 81: {
 82:   TS             ts = (TS)ctx;
 84:   Vec            X0,Xdot,X0_c,Xdot_c;

 87:   TSThetaGetX0AndXdot(ts,fine,&X0,&Xdot);
 88:   TSThetaGetX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
 89:   MatRestrict(restrct,X0,X0_c);
 90:   MatRestrict(restrct,Xdot,Xdot_c);
 91:   VecPointwiseMult(X0_c,rscale,X0_c);
 92:   VecPointwiseMult(Xdot_c,rscale,Xdot_c);
 93:   TSThetaRestoreX0AndXdot(ts,fine,&X0,&Xdot);
 94:   TSThetaRestoreX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
 95:   return(0);
 96: }

100: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm,DM subdm,void *ctx)
101: {

104:   return(0);
105: }

109: static PetscErrorCode DMSubDomainRestrictHook_TSTheta(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
110: {
111:   TS             ts = (TS)ctx;
113:   Vec            X0,Xdot,X0_sub,Xdot_sub;

116:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
117:   TSThetaGetX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);

119:   VecScatterBegin(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
120:   VecScatterEnd(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);

122:   VecScatterBegin(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
123:   VecScatterEnd(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);

125:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
126:   TSThetaRestoreX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
127:   return(0);
128: }

132: static PetscErrorCode TSEvaluateStep_Theta(TS ts,PetscInt order,Vec U,PetscBool *done)
133: {
135:   TS_Theta       *th = (TS_Theta*)ts->data;

138:   if (order == 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"No time-step adaptivity implemented for 1st order theta method; Run with -ts_adapt_type none");
139:   if (order == th->order) {
140:     if (th->endpoint) {
141:       VecCopy(th->X,U);
142:     } else {
143:       PetscReal shift = 1./(th->Theta*ts->time_step);
144:       VecAXPBYPCZ(th->Xdot,-shift,shift,0,U,th->X);
145:       VecAXPY(U,ts->time_step,th->Xdot);
146:     }
147:   } else if (order == th->order-1 && order) {
148:     VecWAXPY(U,ts->time_step,th->Xdot,th->X0);
149:   }
150:   return(0);
151: }

155: static PetscErrorCode TSRollBack_Theta(TS ts)
156: {
157:   TS_Theta       *th = (TS_Theta*)ts->data;

161:   VecCopy(th->X0,ts->vec_sol);
162:   th->status    = TS_STEP_INCOMPLETE;
163:   return(0);
164: }

168: static PetscErrorCode TSStep_Theta(TS ts)
169: {
170:   TS_Theta       *th = (TS_Theta*)ts->data;
171:   PetscInt       its,lits,reject,next_scheme;
172:   PetscReal      next_time_step;
173:   TSAdapt        adapt;
174:   PetscBool      stageok,accept = PETSC_TRUE;

178:   th->status = TS_STEP_INCOMPLETE;
179:   VecCopy(ts->vec_sol,th->X0);
180:   for (reject=0; !ts->reason && th->status != TS_STEP_COMPLETE; ts->reject++) {
181:     PetscReal shift = 1./(th->Theta*ts->time_step);
182:     th->stage_time = ts->ptime + (th->endpoint ? 1. : th->Theta)*ts->time_step;
183:     TSPreStep(ts);
184:     TSPreStage(ts,th->stage_time);

186:     if (th->endpoint) {           /* This formulation assumes linear time-independent mass matrix */
187:       VecZeroEntries(th->Xdot);
188:       if (!th->affine) {VecDuplicate(ts->vec_sol,&th->affine);}
189:       TSComputeIFunction(ts,ts->ptime,ts->vec_sol,th->Xdot,th->affine,PETSC_FALSE);
190:       VecScale(th->affine,(th->Theta-1.)/th->Theta);
191:     }
192:     if (th->extrapolate) {
193:       VecWAXPY(th->X,1./shift,th->Xdot,ts->vec_sol);
194:     } else {
195:       VecCopy(ts->vec_sol,th->X);
196:     }
197:     SNESSolve(ts->snes,th->affine,th->X);
198:     SNESGetIterationNumber(ts->snes,&its);
199:     SNESGetLinearSolveIterations(ts->snes,&lits);
200:     ts->snes_its += its; ts->ksp_its += lits;
201:     TSPostStage(ts,th->stage_time,0,&(th->X));
202:     TSGetAdapt(ts,&adapt);
203:     TSAdaptCheckStage(adapt,ts,th->stage_time,th->X,&stageok);
204:     if (!stageok) {accept = PETSC_FALSE; goto reject_step;}

206:     TSEvaluateStep(ts,th->order,ts->vec_sol,NULL);
207:     th->status = TS_STEP_PENDING;
208:     /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
209:     TSGetAdapt(ts,&adapt);
210:     TSAdaptCandidatesClear(adapt);
211:     TSAdaptCandidateAdd(adapt,NULL,th->order,1,th->ccfl,1.0,PETSC_TRUE);
212:     TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
213:     if (!accept) {           /* Roll back the current step */
214:       ts->ptime += next_time_step; /* This will be undone in rollback */
215:       th->status = TS_STEP_INCOMPLETE;
216:       TSRollBack(ts);
217:       goto reject_step;
218:     }

220:     if (ts->vec_costintegral) {
221:       /* Evolve ts->vec_costintegral to compute integrals */
222:       if (th->endpoint) {
223:         TSAdjointComputeCostIntegrand(ts,ts->ptime,th->X0,ts->vec_costintegrand);
224:         VecAXPY(ts->vec_costintegral,ts->time_step*(1.-th->Theta),ts->vec_costintegrand);
225:       }
226:       TSAdjointComputeCostIntegrand(ts,th->stage_time,th->X,ts->vec_costintegrand);
227:       if (th->endpoint) {
228:         VecAXPY(ts->vec_costintegral,ts->time_step*th->Theta,ts->vec_costintegrand);
229:       }else {
230:         VecAXPY(ts->vec_costintegral,ts->time_step,ts->vec_costintegrand);
231:       }
232:     }

234:     /* ignore next_scheme for now */
235:     ts->ptime    += ts->time_step;
236:     ts->time_step = next_time_step;
237:     ts->steps++;
238:     th->status = TS_STEP_COMPLETE;
239:     break;

241: reject_step:
242:     if (!ts->reason && ++reject > ts->max_reject && ts->max_reject >= 0) {
243:       ts->reason = TS_DIVERGED_STEP_REJECTED;
244:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
245:     }
246:     continue;
247:   }
248:   return(0);
249: }

253: static PetscErrorCode TSAdjointStep_Theta(TS ts)
254: {
255:   TS_Theta            *th = (TS_Theta*)ts->data;
256:   Vec                 *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
257:   PetscInt            nadj;
258:   PetscErrorCode      ierr;
259:   Mat                 J,Jp;
260:   KSP                 ksp;
261:   PetscReal           shift;


265:   th->status = TS_STEP_INCOMPLETE;
266:   SNESGetKSP(ts->snes,&ksp);
267:   TSGetIJacobian(ts,&J,&Jp,NULL,NULL);

269:   /* If endpoint=1, th->ptime and th->X0 will be used; if endpoint=0, th->stage_time and th->X will be used. */
270:   th->stage_time = ts->ptime + (th->endpoint ? ts->time_step : (1.-th->Theta)*ts->time_step); /* time_step is negative*/
271:   th->ptime      = ts->ptime + ts->time_step;

273:   TSPreStep(ts);

275:   /* Build RHS */
276:   if (ts->vec_costintegral) { /* Cost function has an integral term */
277:     if (th->endpoint) {
278:       TSAdjointComputeDRDYFunction(ts,ts->ptime,ts->vec_sol,ts->vecs_drdy);
279:     }else {
280:       TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
281:     }
282:   }
283:   for (nadj=0; nadj<ts->numcost; nadj++) {
284:     VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);
285:     VecScale(VecsSensiTemp[nadj],-1./(th->Theta*ts->time_step));
286:     if (ts->vec_costintegral) {
287:       VecAXPY(VecsSensiTemp[nadj],1.,ts->vecs_drdy[nadj]);
288:     }
289:   }

291:   /* Build LHS */
292:   shift = -1./(th->Theta*ts->time_step);
293:   if (th->endpoint) {
294:     TSComputeIJacobian(ts,ts->ptime,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE);
295:   }else {
296:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
297:   }
298:   KSPSetOperators(ksp,J,Jp);

300:   /* Solve LHS X = RHS */
301:   for (nadj=0; nadj<ts->numcost; nadj++) {
302:     KSPSolveTranspose(ksp,VecsSensiTemp[nadj],VecsDeltaLam[nadj]);
303:   }

305:   /* Update sensitivities, and evaluate integrals if there is any */
306:   if(th->endpoint) { /* two-stage case */
307:     if (th->Theta!=1.) {
308:       shift = -1./((th->Theta-1.)*ts->time_step);
309:       TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,shift,J,Jp,PETSC_FALSE);
310:       if (ts->vec_costintegral) {
311:         TSAdjointComputeDRDYFunction(ts,th->ptime,th->X0,ts->vecs_drdy);
312:         if (!ts->costintegralfwd) {
313:           /* Evolve ts->vec_costintegral to compute integrals */
314:           TSAdjointComputeCostIntegrand(ts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
315:           VecAXPY(ts->vec_costintegral,-ts->time_step*th->Theta,ts->vec_costintegrand);
316:           TSAdjointComputeCostIntegrand(ts,th->ptime,th->X0,ts->vec_costintegrand);
317:           VecAXPY(ts->vec_costintegral,ts->time_step*(th->Theta-1.),ts->vec_costintegrand);
318:         }
319:       }
320:       for (nadj=0; nadj<ts->numcost; nadj++) {
321:         MatMultTranspose(J,VecsDeltaLam[nadj],ts->vecs_sensi[nadj]);
322:         if (ts->vec_costintegral) {
323:           VecAXPY(ts->vecs_sensi[nadj],-1.,ts->vecs_drdy[nadj]);
324:         }
325:         VecScale(ts->vecs_sensi[nadj],1./shift);
326:       }
327:     }else { /* backward Euler */
328:       shift = 0.0;
329:       TSComputeIJacobian(ts,ts->ptime,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
330:       for (nadj=0; nadj<ts->numcost; nadj++) {
331:         MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
332:         VecAXPY(ts->vecs_sensi[nadj],ts->time_step,VecsSensiTemp[nadj]);
333:         if (ts->vec_costintegral) {
334:           VecAXPY(ts->vecs_sensi[nadj],-ts->time_step,ts->vecs_drdy[nadj]);
335:           if (!ts->costintegralfwd) {
336:           /* Evolve ts->vec_costintegral to compute integrals */
337:             TSAdjointComputeCostIntegrand(ts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
338:             VecAXPY(ts->vec_costintegral,-ts->time_step*th->Theta,ts->vec_costintegrand);
339:           }
340:         }
341:       }
342:     }

344:     if (ts->vecs_sensip) { /* sensitivities wrt parameters */
345:       TSAdjointComputeRHSJacobian(ts,ts->ptime,ts->vec_sol,ts->Jacp);
346:       for (nadj=0; nadj<ts->numcost; nadj++) {
347:         MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
348:         VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*th->Theta,VecsDeltaMu[nadj]);
349:       }
350:       if (th->Theta!=1.) {
351:         TSAdjointComputeRHSJacobian(ts,th->ptime,th->X0,ts->Jacp);
352:         for (nadj=0; nadj<ts->numcost; nadj++) {
353:           MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
354:           VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*(1.-th->Theta),VecsDeltaMu[nadj]);
355:         }
356:       }
357:       if (ts->vec_costintegral) {
358:         TSAdjointComputeDRDPFunction(ts,ts->ptime,ts->vec_sol,ts->vecs_drdp);
359:         for (nadj=0; nadj<ts->numcost; nadj++) {
360:           VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*th->Theta,ts->vecs_drdp[nadj]);
361:         }
362:         if (th->Theta!=1.) {
363:           TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
364:           for (nadj=0; nadj<ts->numcost; nadj++) {
365:             VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*(1.-th->Theta),ts->vecs_drdp[nadj]);
366:           }
367:         }
368:       }
369:     }
370:   }else { /* one-stage case */
371:     shift = 0.0;
372:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
373:     if (ts->vec_costintegral) {
374:       TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
375:       if (!ts->costintegralfwd) {
376:         /* Evolve ts->vec_costintegral to compute integrals */
377:         TSAdjointComputeCostIntegrand(ts,th->stage_time,th->X,ts->vec_costintegrand);
378:         VecAXPY(ts->vec_costintegral,-ts->time_step,ts->vec_costintegrand);
379:       }
380:     }
381:     for (nadj=0; nadj<ts->numcost; nadj++) {
382:       MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
383:       VecAXPY(ts->vecs_sensi[nadj],ts->time_step,VecsSensiTemp[nadj]);
384:       if (ts->vec_costintegral) {
385:         VecAXPY(ts->vecs_sensi[nadj],-ts->time_step,ts->vecs_drdy[nadj]);
386:       }
387:     }
388:     if (ts->vecs_sensip) {
389:       TSAdjointComputeRHSJacobian(ts,th->stage_time,th->X,ts->Jacp);
390:       for (nadj=0; nadj<ts->numcost; nadj++) {
391:         MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
392:         VecAXPY(ts->vecs_sensip[nadj],-ts->time_step,VecsDeltaMu[nadj]);
393:       }
394:       if (ts->vec_costintegral) {
395:         TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
396:         for (nadj=0; nadj<ts->numcost; nadj++) {
397:           VecAXPY(ts->vecs_sensip[nadj],-ts->time_step,ts->vecs_drdp[nadj]);
398:         }
399:       }
400:     }
401:   }

403:   ts->ptime += ts->time_step;
404:   ts->steps++;
405:   th->status = TS_STEP_COMPLETE;
406:   return(0);
407: }

411: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
412: {
413:   TS_Theta       *th   = (TS_Theta*)ts->data;
414:   PetscReal      alpha = t - ts->ptime;

418:   VecCopy(ts->vec_sol,th->X);
419:   if (th->endpoint) alpha *= th->Theta;
420:   VecWAXPY(X,alpha,th->Xdot,th->X);
421:   return(0);
422: }

424: /*------------------------------------------------------------*/
427: static PetscErrorCode TSReset_Theta(TS ts)
428: {
429:   TS_Theta       *th = (TS_Theta*)ts->data;

433:   VecDestroy(&th->X);
434:   VecDestroy(&th->Xdot);
435:   VecDestroy(&th->X0);
436:   VecDestroy(&th->affine);
437:   VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
438:   VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
439:   VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);
440:   return(0);
441: }

445: static PetscErrorCode TSDestroy_Theta(TS ts)
446: {

450:   TSReset_Theta(ts);
451:   PetscFree(ts->data);
452:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
453:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
454:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
455:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
456:   return(0);
457: }

459: /*
460:   This defines the nonlinear equation that is to be solved with SNES
461:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
462: */
465: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
466: {
467:   TS_Theta       *th = (TS_Theta*)ts->data;
469:   Vec            X0,Xdot;
470:   DM             dm,dmsave;
471:   PetscReal      shift = 1./(th->Theta*ts->time_step);

474:   SNESGetDM(snes,&dm);
475:   /* When using the endpoint variant, this is actually 1/Theta * Xdot */
476:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
477:   VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);

479:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
480:   dmsave = ts->dm;
481:   ts->dm = dm;
482:   TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
483:   ts->dm = dmsave;
484:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
485:   return(0);
486: }

490: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
491: {
492:   TS_Theta       *th = (TS_Theta*)ts->data;
494:   Vec            Xdot;
495:   DM             dm,dmsave;
496:   PetscReal      shift = 1./(th->Theta*ts->time_step);

499:   SNESGetDM(snes,&dm);

501:   /* th->Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
502:   TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);

504:   dmsave = ts->dm;
505:   ts->dm = dm;
506:   TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
507:   ts->dm = dmsave;
508:   TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
509:   return(0);
510: }

514: static PetscErrorCode TSSetUp_Theta(TS ts)
515: {
516:   TS_Theta       *th = (TS_Theta*)ts->data;
518:   SNES           snes;
519:   TSAdapt        adapt;
520:   DM             dm;

523:   if (!th->X) {
524:     VecDuplicate(ts->vec_sol,&th->X);
525:   }
526:   if (!th->Xdot) {
527:     VecDuplicate(ts->vec_sol,&th->Xdot);
528:   }
529:   if (!th->X0) {
530:     VecDuplicate(ts->vec_sol,&th->X0);
531:   }
532:   TSGetSNES(ts,&snes);
533:   TSGetDM(ts,&dm);
534:   if (dm) {
535:     DMCoarsenHookAdd(dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
536:     DMSubDomainHookAdd(dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
537:   }
538:   if (th->Theta == 0.5 && th->endpoint) th->order = 2;
539:   else th->order = 1;

541:   TSGetAdapt(ts,&adapt);
542:   if (!th->adapt) {
543:     TSAdaptSetType(adapt,TSADAPTNONE);
544:   }
545:   return(0);
546: }

550: static PetscErrorCode TSSetUp_BEuler(TS ts)
551: {
552:   TS_Theta       *th = (TS_Theta*)ts->data;

556:   if (th->Theta != 1.0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change the default value (1) of theta when using backward Euler\n");
557:   TSSetUp_Theta(ts);
558:   return(0);
559: }

563: static PetscErrorCode TSSetUp_CN(TS ts)
564: {
565:   TS_Theta       *th = (TS_Theta*)ts->data;

569:   if (th->Theta != 0.5) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change the default value (0.5) of theta when using Crank-Nicolson\n");
570:   if (!th->endpoint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change to the midpoint form of the Theta methods when using Crank-Nicolson\n");
571:   TSSetUp_Theta(ts);
572:   return(0);
573: }
574: /*------------------------------------------------------------*/

578: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
579: {
580:   TS_Theta       *th = (TS_Theta*)ts->data;

584:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
585:   if(ts->vecs_sensip) {
586:     VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
587:   }
588:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
589:   return(0);
590: }
591: /*------------------------------------------------------------*/

595: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
596: {
597:   TS_Theta       *th = (TS_Theta*)ts->data;

601:   PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
602:   {
603:     PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
604:     PetscOptionsBool("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
605:     PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
606:     PetscOptionsBool("-ts_theta_adapt","Use time-step adaptivity with the Theta method","",th->adapt,&th->adapt,NULL);
607:   }
608:   PetscOptionsTail();
609:   return(0);
610: }

614: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
615: {
616:   TS_Theta       *th = (TS_Theta*)ts->data;
617:   PetscBool      iascii;

621:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
622:   if (iascii) {
623:     PetscViewerASCIIPrintf(viewer,"  Theta=%g\n",(double)th->Theta);
624:     PetscViewerASCIIPrintf(viewer,"  Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
625:   }
626:   if (ts->snes) {SNESView(ts->snes,viewer);}
627:   return(0);
628: }

632: PetscErrorCode  TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
633: {
634:   TS_Theta *th = (TS_Theta*)ts->data;

637:   *theta = th->Theta;
638:   return(0);
639: }

643: PetscErrorCode  TSThetaSetTheta_Theta(TS ts,PetscReal theta)
644: {
645:   TS_Theta *th = (TS_Theta*)ts->data;

648:   if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
649:   th->Theta = theta;
650:   return(0);
651: }

655: PetscErrorCode  TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
656: {
657:   TS_Theta *th = (TS_Theta*)ts->data;

660:   *endpoint = th->endpoint;
661:   return(0);
662: }

666: PetscErrorCode  TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
667: {
668:   TS_Theta *th = (TS_Theta*)ts->data;

671:   th->endpoint = flg;
672:   return(0);
673: }

675: #if defined(PETSC_HAVE_COMPLEX)
678: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
679: {
680:   PetscComplex z   = xr + xi*PETSC_i,f;
681:   TS_Theta     *th = (TS_Theta*)ts->data;
682:   const PetscReal one = 1.0;

685:   f   = (one + (one - th->Theta)*z)/(one - th->Theta*z);
686:   *yr = PetscRealPartComplex(f);
687:   *yi = PetscImaginaryPartComplex(f);
688:   return(0);
689: }
690: #endif

694: static PetscErrorCode  TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
695: {
696:   TS_Theta     *th = (TS_Theta*)ts->data;

699:   *ns = 1;
700:   if(Y) {
701:     *Y  = (th->endpoint)?&(th->X0):&(th->X);
702:   }
703:   return(0);
704: }

706: /* ------------------------------------------------------------ */
707: /*MC
708:       TSTHETA - DAE solver using the implicit Theta method

710:    Level: beginner

712:    Options Database:
713: +      -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
714: .      -ts_theta_extrapolate <flg> - Extrapolate stage solution from previous solution (sometimes unstable)
715: .      -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of the Theta method
716: -     -ts_theta_adapt <flg> - Use time-step adaptivity with the Theta method

718:    Notes:
719: $  -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
720: $  -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
721: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)

723:    This method can be applied to DAE.

725:    This method is cast as a 1-stage implicit Runge-Kutta method.

727: .vb
728:   Theta | Theta
729:   -------------
730:         |  1
731: .ve

733:    For the default Theta=0.5, this is also known as the implicit midpoint rule.

735:    When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:

737: .vb
738:   0 | 0         0
739:   1 | 1-Theta   Theta
740:   -------------------
741:     | 1-Theta   Theta
742: .ve

744:    For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).

746:    To apply a diagonally implicit RK method to DAE, the stage formula

748: $  Y_i = X + h sum_j a_ij Y'_j

750:    is interpreted as a formula for Y'_i in terms of Y_i and known values (Y'_j, j<i)

752: .seealso:  TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()

754: M*/
757: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
758: {
759:   TS_Theta       *th;

763:   ts->ops->reset           = TSReset_Theta;
764:   ts->ops->destroy         = TSDestroy_Theta;
765:   ts->ops->view            = TSView_Theta;
766:   ts->ops->setup           = TSSetUp_Theta;
767:   ts->ops->adjointsetup    = TSAdjointSetUp_Theta;
768:   ts->ops->step            = TSStep_Theta;
769:   ts->ops->interpolate     = TSInterpolate_Theta;
770:   ts->ops->evaluatestep    = TSEvaluateStep_Theta;
771:   ts->ops->rollback        = TSRollBack_Theta;
772:   ts->ops->setfromoptions  = TSSetFromOptions_Theta;
773:   ts->ops->snesfunction    = SNESTSFormFunction_Theta;
774:   ts->ops->snesjacobian    = SNESTSFormJacobian_Theta;
775: #if defined(PETSC_HAVE_COMPLEX)
776:   ts->ops->linearstability = TSComputeLinearStability_Theta;
777: #endif
778:   ts->ops->getstages       = TSGetStages_Theta;
779:   ts->ops->adjointstep     = TSAdjointStep_Theta;

781:   PetscNewLog(ts,&th);
782:   ts->data = (void*)th;

784:   th->extrapolate = PETSC_FALSE;
785:   th->Theta       = 0.5;
786:   th->ccfl        = 1.0;
787:   th->adapt       = PETSC_FALSE;
788:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
789:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
790:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
791:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
792:   return(0);
793: }

797: /*@
798:   TSThetaGetTheta - Get the abscissa of the stage in (0,1].

800:   Not Collective

802:   Input Parameter:
803: .  ts - timestepping context

805:   Output Parameter:
806: .  theta - stage abscissa

808:   Note:
809:   Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.

811:   Level: Advanced

813: .seealso: TSThetaSetTheta()
814: @*/
815: PetscErrorCode  TSThetaGetTheta(TS ts,PetscReal *theta)
816: {

822:   PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
823:   return(0);
824: }

828: /*@
829:   TSThetaSetTheta - Set the abscissa of the stage in (0,1].

831:   Not Collective

833:   Input Parameter:
834: +  ts - timestepping context
835: -  theta - stage abscissa

837:   Options Database:
838: .  -ts_theta_theta <theta>

840:   Level: Intermediate

842: .seealso: TSThetaGetTheta()
843: @*/
844: PetscErrorCode  TSThetaSetTheta(TS ts,PetscReal theta)
845: {

850:   PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
851:   return(0);
852: }

856: /*@
857:   TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).

859:   Not Collective

861:   Input Parameter:
862: .  ts - timestepping context

864:   Output Parameter:
865: .  endpoint - PETSC_TRUE when using the endpoint variant

867:   Level: Advanced

869: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
870: @*/
871: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
872: {

878:   PetscTryMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
879:   return(0);
880: }

884: /*@
885:   TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).

887:   Not Collective

889:   Input Parameter:
890: +  ts - timestepping context
891: -  flg - PETSC_TRUE to use the endpoint variant

893:   Options Database:
894: .  -ts_theta_endpoint <flg>

896:   Level: Intermediate

898: .seealso: TSTHETA, TSCN
899: @*/
900: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
901: {

906:   PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
907:   return(0);
908: }

910: /*
911:  * TSBEULER and TSCN are straightforward specializations of TSTHETA.
912:  * The creation functions for these specializations are below.
913:  */

917: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
918: {

922:   SNESView(ts->snes,viewer);
923:   return(0);
924: }

926: /*MC
927:       TSBEULER - ODE solver using the implicit backward Euler method

929:   Level: beginner

931:   Notes:
932:   TSBEULER is equivalent to TSTHETA with Theta=1.0

934: $  -ts_type theta -ts_theta_theta 1.

936: .seealso:  TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA

938: M*/
941: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
942: {

946:   TSCreate_Theta(ts);
947:   TSThetaSetTheta(ts,1.0);
948:   ts->ops->setup = TSSetUp_BEuler;
949:   ts->ops->view = TSView_BEuler;
950:   return(0);
951: }

955: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
956: {

960:   SNESView(ts->snes,viewer);
961:   return(0);
962: }

964: /*MC
965:       TSCN - ODE solver using the implicit Crank-Nicolson method.

967:   Level: beginner

969:   Notes:
970:   TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.

972: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

974: .seealso:  TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA

976: M*/
979: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
980: {

984:   TSCreate_Theta(ts);
985:   TSThetaSetTheta(ts,0.5);
986:   TSThetaSetEndpoint(ts,PETSC_TRUE);
987:   ts->ops->setup = TSSetUp_CN;
988:   ts->ops->view = TSView_CN;
989:   return(0);
990: }