Actual source code: theta.c
petsc-3.6.4 2016-04-12
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: }
52: static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
53: {
57: if (X0) {
58: if (dm && dm != ts->dm) {
59: DMRestoreNamedGlobalVector(dm,"TSTheta_X0",X0);
60: }
61: }
62: if (Xdot) {
63: if (dm && dm != ts->dm) {
64: DMRestoreNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
65: }
66: }
67: return(0);
68: }
72: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine,DM coarse,void *ctx)
73: {
76: return(0);
77: }
81: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
82: {
83: TS ts = (TS)ctx;
85: Vec X0,Xdot,X0_c,Xdot_c;
88: TSThetaGetX0AndXdot(ts,fine,&X0,&Xdot);
89: TSThetaGetX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
90: MatRestrict(restrct,X0,X0_c);
91: MatRestrict(restrct,Xdot,Xdot_c);
92: VecPointwiseMult(X0_c,rscale,X0_c);
93: VecPointwiseMult(Xdot_c,rscale,Xdot_c);
94: TSThetaRestoreX0AndXdot(ts,fine,&X0,&Xdot);
95: TSThetaRestoreX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
96: return(0);
97: }
101: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm,DM subdm,void *ctx)
102: {
105: return(0);
106: }
110: static PetscErrorCode DMSubDomainRestrictHook_TSTheta(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
111: {
112: TS ts = (TS)ctx;
114: Vec X0,Xdot,X0_sub,Xdot_sub;
117: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
118: TSThetaGetX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
120: VecScatterBegin(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
121: VecScatterEnd(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
123: VecScatterBegin(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
124: VecScatterEnd(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
126: TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
127: TSThetaRestoreX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
128: return(0);
129: }
133: static PetscErrorCode TSEvaluateStep_Theta(TS ts,PetscInt order,Vec U,PetscBool *done)
134: {
136: TS_Theta *th = (TS_Theta*)ts->data;
139: 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");
140: if (order == th->order) {
141: if (th->endpoint) {
142: VecCopy(th->X,U);
143: } else {
144: PetscReal shift = 1./(th->Theta*ts->time_step);
145: VecAXPBYPCZ(th->Xdot,-shift,shift,0,U,th->X);
146: VecAXPY(U,ts->time_step,th->Xdot);
147: }
148: } else if (order == th->order-1 && order) {
149: VecWAXPY(U,ts->time_step,th->Xdot,th->X0);
150: }
151: return(0);
152: }
156: static PetscErrorCode TSRollBack_Theta(TS ts)
157: {
158: TS_Theta *th = (TS_Theta*)ts->data;
162: VecCopy(th->X0,ts->vec_sol);
163: th->status = TS_STEP_INCOMPLETE;
164: return(0);
165: }
169: static PetscErrorCode TSStep_Theta(TS ts)
170: {
171: TS_Theta *th = (TS_Theta*)ts->data;
172: PetscInt its,lits,reject,next_scheme;
173: PetscReal next_time_step;
174: TSAdapt adapt;
175: PetscBool stageok,accept = PETSC_TRUE;
179: th->status = TS_STEP_INCOMPLETE;
180: VecCopy(ts->vec_sol,th->X0);
181: for (reject=0; !ts->reason && th->status != TS_STEP_COMPLETE; ts->reject++) {
182: PetscReal shift = 1./(th->Theta*ts->time_step);
183: th->stage_time = ts->ptime + (th->endpoint ? 1. : th->Theta)*ts->time_step;
184: TSPreStep(ts);
185: TSPreStage(ts,th->stage_time);
187: if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */
188: VecZeroEntries(th->Xdot);
189: if (!th->affine) {VecDuplicate(ts->vec_sol,&th->affine);}
190: TSComputeIFunction(ts,ts->ptime,ts->vec_sol,th->Xdot,th->affine,PETSC_FALSE);
191: VecScale(th->affine,(th->Theta-1.)/th->Theta);
192: }
193: if (th->extrapolate) {
194: VecWAXPY(th->X,1./shift,th->Xdot,ts->vec_sol);
195: } else {
196: VecCopy(ts->vec_sol,th->X);
197: }
198: SNESSolve(ts->snes,th->affine,th->X);
199: SNESGetIterationNumber(ts->snes,&its);
200: SNESGetLinearSolveIterations(ts->snes,&lits);
201: ts->snes_its += its; ts->ksp_its += lits;
202: TSPostStage(ts,th->stage_time,0,&(th->X));
203: TSGetAdapt(ts,&adapt);
204: TSAdaptCheckStage(adapt,ts,&stageok);
205: if (!stageok) {accept = PETSC_FALSE; goto reject_step;}
207: TSEvaluateStep(ts,th->order,ts->vec_sol,NULL);
208: th->status = TS_STEP_PENDING;
209: /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */
210: TSGetAdapt(ts,&adapt);
211: TSAdaptCandidatesClear(adapt);
212: TSAdaptCandidateAdd(adapt,NULL,th->order,1,th->ccfl,1.0,PETSC_TRUE);
213: TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);
214: if (!accept) { /* Roll back the current step */
215: ts->ptime += next_time_step; /* This will be undone in rollback */
216: th->status = TS_STEP_INCOMPLETE;
217: TSRollBack(ts);
218: goto reject_step;
219: }
221: if (ts->vec_costintegral) {
222: /* Evolve ts->vec_costintegral to compute integrals */
223: if (th->endpoint) {
224: TSAdjointComputeCostIntegrand(ts,ts->ptime,th->X0,ts->vec_costintegrand);
225: VecAXPY(ts->vec_costintegral,ts->time_step*(1.-th->Theta),ts->vec_costintegrand);
226: }
227: TSAdjointComputeCostIntegrand(ts,th->stage_time,th->X,ts->vec_costintegrand);
228: if (th->endpoint) {
229: VecAXPY(ts->vec_costintegral,ts->time_step*th->Theta,ts->vec_costintegrand);
230: }else {
231: VecAXPY(ts->vec_costintegral,ts->time_step,ts->vec_costintegrand);
232: }
233: }
235: /* ignore next_scheme for now */
236: ts->ptime += ts->time_step;
237: ts->time_step = next_time_step;
238: ts->steps++;
239: th->status = TS_STEP_COMPLETE;
240: break;
242: reject_step:
243: if (!ts->reason && ++reject > ts->max_reject && ts->max_reject >= 0) {
244: ts->reason = TS_DIVERGED_STEP_REJECTED;
245: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,reject);
246: }
247: continue;
248: }
249: return(0);
250: }
254: static PetscErrorCode TSAdjointStep_Theta(TS ts)
255: {
256: TS_Theta *th = (TS_Theta*)ts->data;
257: Vec *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
258: PetscInt nadj;
259: PetscErrorCode ierr;
260: Mat J,Jp;
261: KSP ksp;
262: PetscReal shift;
265:
266: th->status = TS_STEP_INCOMPLETE;
267: SNESGetKSP(ts->snes,&ksp);
268: TSGetIJacobian(ts,&J,&Jp,NULL,NULL);
270: /* If endpoint=1, th->ptime and th->X0 will be used; if endpoint=0, th->stage_time and th->X will be used. */
271: th->stage_time = ts->ptime + (th->endpoint ? ts->time_step : (1.-th->Theta)*ts->time_step); /* time_step is negative*/
272: th->ptime = ts->ptime + ts->time_step;
274: TSPreStep(ts);
276: /* Build RHS */
277: if (ts->vec_costintegral) { /* Cost function has an integral term */
278: if (th->endpoint) {
279: TSAdjointComputeDRDYFunction(ts,ts->ptime,ts->vec_sol,ts->vecs_drdy);
280: }else {
281: TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
282: }
283: }
284: for (nadj=0; nadj<ts->numcost; nadj++) {
285: VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);
286: VecScale(VecsSensiTemp[nadj],-1./(th->Theta*ts->time_step));
287: if (ts->vec_costintegral) {
288: VecAXPY(VecsSensiTemp[nadj],1.,ts->vecs_drdy[nadj]);
289: }
290: }
291:
292: /* Build LHS */
293: shift = -1./(th->Theta*ts->time_step);
294: if (th->endpoint) {
295: TSComputeIJacobian(ts,ts->ptime,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE);
296: }else {
297: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
298: }
299: KSPSetOperators(ksp,J,Jp);
301: /* Solve LHS X = RHS */
302: for (nadj=0; nadj<ts->numcost; nadj++) {
303: KSPSolveTranspose(ksp,VecsSensiTemp[nadj],VecsDeltaLam[nadj]);
304: }
306: /* Update sensitivities, and evaluate integrals if there is any */
307: if(th->endpoint) { /* two-stage case */
308: if (th->Theta!=1.) {
309: shift = -1./((th->Theta-1.)*ts->time_step);
310: TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,shift,J,Jp,PETSC_FALSE);
311: if (ts->vec_costintegral) {
312: TSAdjointComputeDRDYFunction(ts,th->ptime,th->X0,ts->vecs_drdy);
313: if (!ts->costintegralfwd) {
314: /* Evolve ts->vec_costintegral to compute integrals */
315: TSAdjointComputeCostIntegrand(ts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
316: VecAXPY(ts->vec_costintegral,-ts->time_step*th->Theta,ts->vec_costintegrand);
317: TSAdjointComputeCostIntegrand(ts,th->ptime,th->X0,ts->vec_costintegrand);
318: VecAXPY(ts->vec_costintegral,ts->time_step*(th->Theta-1.),ts->vec_costintegrand);
319: }
320: }
321: for (nadj=0; nadj<ts->numcost; nadj++) {
322: MatMultTranspose(J,VecsDeltaLam[nadj],ts->vecs_sensi[nadj]);
323: if (ts->vec_costintegral) {
324: VecAXPY(ts->vecs_sensi[nadj],-1.,ts->vecs_drdy[nadj]);
325: }
326: VecScale(ts->vecs_sensi[nadj],1./shift);
327: }
328: }else { /* backward Euler */
329: shift = 0.0;
330: TSComputeIJacobian(ts,ts->ptime,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
331: for (nadj=0; nadj<ts->numcost; nadj++) {
332: MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
333: VecAXPY(ts->vecs_sensi[nadj],ts->time_step,VecsSensiTemp[nadj]);
334: if (ts->vec_costintegral) {
335: VecAXPY(ts->vecs_sensi[nadj],-ts->time_step,ts->vecs_drdy[nadj]);
336: if (!ts->costintegralfwd) {
337: /* Evolve ts->vec_costintegral to compute integrals */
338: TSAdjointComputeCostIntegrand(ts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
339: VecAXPY(ts->vec_costintegral,-ts->time_step*th->Theta,ts->vec_costintegrand);
340: }
341: }
342: }
343: }
345: if (ts->vecs_sensip) { /* sensitivities wrt parameters */
346: TSAdjointComputeRHSJacobian(ts,ts->ptime,ts->vec_sol,ts->Jacp);
347: for (nadj=0; nadj<ts->numcost; nadj++) {
348: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
349: VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*th->Theta,VecsDeltaMu[nadj]);
350: }
351: if (th->Theta!=1.) {
352: TSAdjointComputeRHSJacobian(ts,th->ptime,th->X0,ts->Jacp);
353: for (nadj=0; nadj<ts->numcost; nadj++) {
354: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
355: VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*(1.-th->Theta),VecsDeltaMu[nadj]);
356: }
357: }
358: if (ts->vec_costintegral) {
359: TSAdjointComputeDRDPFunction(ts,ts->ptime,ts->vec_sol,ts->vecs_drdp);
360: for (nadj=0; nadj<ts->numcost; nadj++) {
361: VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*th->Theta,ts->vecs_drdp[nadj]);
362: }
363: if (th->Theta!=1.) {
364: TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
365: for (nadj=0; nadj<ts->numcost; nadj++) {
366: VecAXPY(ts->vecs_sensip[nadj],-ts->time_step*(1.-th->Theta),ts->vecs_drdp[nadj]);
367: }
368: }
369: }
370: }
371: }else { /* one-stage case */
372: shift = 0.0;
373: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
374: if (ts->vec_costintegral) {
375: TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
376: if (!ts->costintegralfwd) {
377: /* Evolve ts->vec_costintegral to compute integrals */
378: TSAdjointComputeCostIntegrand(ts,th->stage_time,th->X,ts->vec_costintegrand);
379: VecAXPY(ts->vec_costintegral,-ts->time_step,ts->vec_costintegrand);
380: }
381: }
382: for (nadj=0; nadj<ts->numcost; nadj++) {
383: MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
384: VecAXPY(ts->vecs_sensi[nadj],ts->time_step,VecsSensiTemp[nadj]);
385: if (ts->vec_costintegral) {
386: VecAXPY(ts->vecs_sensi[nadj],-ts->time_step,ts->vecs_drdy[nadj]);
387: }
388: }
389: if (ts->vecs_sensip) {
390: TSAdjointComputeRHSJacobian(ts,th->stage_time,th->X,ts->Jacp);
391: for (nadj=0; nadj<ts->numcost; nadj++) {
392: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
393: VecAXPY(ts->vecs_sensip[nadj],-ts->time_step,VecsDeltaMu[nadj]);
394: }
395: if (ts->vec_costintegral) {
396: TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
397: for (nadj=0; nadj<ts->numcost; nadj++) {
398: VecAXPY(ts->vecs_sensip[nadj],-ts->time_step,ts->vecs_drdp[nadj]);
399: }
400: }
401: }
402: }
403:
404: ts->ptime += ts->time_step;
405: ts->steps++;
406: th->status = TS_STEP_COMPLETE;
407: return(0);
408: }
412: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
413: {
414: TS_Theta *th = (TS_Theta*)ts->data;
415: PetscReal alpha = t - ts->ptime;
419: VecCopy(ts->vec_sol,th->X);
420: if (th->endpoint) alpha *= th->Theta;
421: VecWAXPY(X,alpha,th->Xdot,th->X);
422: return(0);
423: }
425: /*------------------------------------------------------------*/
428: static PetscErrorCode TSReset_Theta(TS ts)
429: {
430: TS_Theta *th = (TS_Theta*)ts->data;
434: VecDestroy(&th->X);
435: VecDestroy(&th->Xdot);
436: VecDestroy(&th->X0);
437: VecDestroy(&th->affine);
438: VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
439: VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
440: VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);
441: return(0);
442: }
446: static PetscErrorCode TSDestroy_Theta(TS ts)
447: {
451: TSReset_Theta(ts);
452: PetscFree(ts->data);
453: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
454: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
455: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
456: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
457: return(0);
458: }
460: /*
461: This defines the nonlinear equation that is to be solved with SNES
462: G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
463: */
466: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
467: {
468: TS_Theta *th = (TS_Theta*)ts->data;
470: Vec X0,Xdot;
471: DM dm,dmsave;
472: PetscReal shift = 1./(th->Theta*ts->time_step);
475: SNESGetDM(snes,&dm);
476: /* When using the endpoint variant, this is actually 1/Theta * Xdot */
477: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
478: VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);
480: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
481: dmsave = ts->dm;
482: ts->dm = dm;
483: TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
484: ts->dm = dmsave;
485: TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
486: return(0);
487: }
491: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
492: {
493: TS_Theta *th = (TS_Theta*)ts->data;
495: Vec Xdot;
496: DM dm,dmsave;
497: PetscReal shift = 1./(th->Theta*ts->time_step);
500: SNESGetDM(snes,&dm);
502: /* th->Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
503: TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);
505: dmsave = ts->dm;
506: ts->dm = dm;
507: TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
508: ts->dm = dmsave;
509: TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
510: return(0);
511: }
515: static PetscErrorCode TSSetUp_Theta(TS ts)
516: {
517: TS_Theta *th = (TS_Theta*)ts->data;
519: SNES snes;
520: TSAdapt adapt;
521: DM dm;
524: if (!th->X) {
525: VecDuplicate(ts->vec_sol,&th->X);
526: }
527: if (!th->Xdot) {
528: VecDuplicate(ts->vec_sol,&th->Xdot);
529: }
530: if (!th->X0) {
531: VecDuplicate(ts->vec_sol,&th->X0);
532: }
533: TSGetSNES(ts,&snes);
534: TSGetDM(ts,&dm);
535: if (dm) {
536: DMCoarsenHookAdd(dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
537: DMSubDomainHookAdd(dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
538: }
539: if (th->Theta == 0.5 && th->endpoint) th->order = 2;
540: else th->order = 1;
542: TSGetAdapt(ts,&adapt);
543: if (!th->adapt) {
544: TSAdaptSetType(adapt,TSADAPTNONE);
545: }
546: return(0);
547: }
551: static PetscErrorCode TSSetUp_BEuler(TS ts)
552: {
553: TS_Theta *th = (TS_Theta*)ts->data;
557: 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");
558: TSSetUp_Theta(ts);
559: return(0);
560: }
564: static PetscErrorCode TSSetUp_CN(TS ts)
565: {
566: TS_Theta *th = (TS_Theta*)ts->data;
570: 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");
571: 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");
572: TSSetUp_Theta(ts);
573: return(0);
574: }
575: /*------------------------------------------------------------*/
579: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
580: {
581: TS_Theta *th = (TS_Theta*)ts->data;
585: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
586: if(ts->vecs_sensip) {
587: VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
588: }
589: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
590: return(0);
591: }
592: /*------------------------------------------------------------*/
596: static PetscErrorCode TSSetFromOptions_Theta(PetscOptions *PetscOptionsObject,TS ts)
597: {
598: TS_Theta *th = (TS_Theta*)ts->data;
602: PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
603: {
604: PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
605: PetscOptionsBool("-ts_theta_extrapolate","Extrapolate stage solution from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
606: PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
607: PetscOptionsBool("-ts_theta_adapt","Use time-step adaptivity with the Theta method","",th->adapt,&th->adapt,NULL);
608: }
609: PetscOptionsTail();
610: return(0);
611: }
615: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
616: {
617: TS_Theta *th = (TS_Theta*)ts->data;
618: PetscBool iascii;
622: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
623: if (iascii) {
624: PetscViewerASCIIPrintf(viewer," Theta=%g\n",(double)th->Theta);
625: PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
626: }
627: if (ts->snes) {SNESView(ts->snes,viewer);}
628: return(0);
629: }
633: PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
634: {
635: TS_Theta *th = (TS_Theta*)ts->data;
638: *theta = th->Theta;
639: return(0);
640: }
644: PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
645: {
646: TS_Theta *th = (TS_Theta*)ts->data;
649: if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
650: th->Theta = theta;
651: return(0);
652: }
656: PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
657: {
658: TS_Theta *th = (TS_Theta*)ts->data;
661: *endpoint = th->endpoint;
662: return(0);
663: }
667: PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
668: {
669: TS_Theta *th = (TS_Theta*)ts->data;
672: th->endpoint = flg;
673: return(0);
674: }
676: #if defined(PETSC_HAVE_COMPLEX)
679: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
680: {
681: PetscComplex z = xr + xi*PETSC_i,f;
682: TS_Theta *th = (TS_Theta*)ts->data;
683: const PetscReal one = 1.0;
686: f = (one + (one - th->Theta)*z)/(one - th->Theta*z);
687: *yr = PetscRealPartComplex(f);
688: *yi = PetscImaginaryPartComplex(f);
689: return(0);
690: }
691: #endif
695: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
696: {
697: TS_Theta *th = (TS_Theta*)ts->data;
700: *ns = 1;
701: if(Y) {
702: *Y = (th->endpoint)?&(th->X0):&(th->X);
703: }
704: return(0);
705: }
707: /* ------------------------------------------------------------ */
708: /*MC
709: TSTHETA - DAE solver using the implicit Theta method
711: Level: beginner
713: Options Database:
714: -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
715: -ts_theta_extrapolate <flg> Extrapolate stage solution from previous solution (sometimes unstable)
716: -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of 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)
725: This method can be applied to DAE.
727: This method is cast as a 1-stage implicit Runge-Kutta method.
729: .vb
730: Theta | Theta
731: -------------
732: | 1
733: .ve
735: For the default Theta=0.5, this is also known as the implicit midpoint rule.
737: When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:
739: .vb
740: 0 | 0 0
741: 1 | 1-Theta Theta
742: -------------------
743: | 1-Theta Theta
744: .ve
746: For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).
748: To apply a diagonally implicit RK method to DAE, the stage formula
750: $ Y_i = X + h sum_j a_ij Y'_j
752: is interpreted as a formula for Y'_i in terms of Y_i and known values (Y'_j, j<i)
754: .seealso: TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()
756: M*/
759: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
760: {
761: TS_Theta *th;
765: ts->ops->reset = TSReset_Theta;
766: ts->ops->destroy = TSDestroy_Theta;
767: ts->ops->view = TSView_Theta;
768: ts->ops->setup = TSSetUp_Theta;
769: ts->ops->adjointsetup = TSAdjointSetUp_Theta;
770: ts->ops->step = TSStep_Theta;
771: ts->ops->interpolate = TSInterpolate_Theta;
772: ts->ops->evaluatestep = TSEvaluateStep_Theta;
773: ts->ops->rollback = TSRollBack_Theta;
774: ts->ops->setfromoptions = TSSetFromOptions_Theta;
775: ts->ops->snesfunction = SNESTSFormFunction_Theta;
776: ts->ops->snesjacobian = SNESTSFormJacobian_Theta;
777: #if defined(PETSC_HAVE_COMPLEX)
778: ts->ops->linearstability = TSComputeLinearStability_Theta;
779: #endif
780: ts->ops->getstages = TSGetStages_Theta;
781: ts->ops->adjointstep = TSAdjointStep_Theta;
783: PetscNewLog(ts,&th);
784: ts->data = (void*)th;
786: th->extrapolate = PETSC_FALSE;
787: th->Theta = 0.5;
788: th->ccfl = 1.0;
789: th->adapt = PETSC_FALSE;
790: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
791: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
792: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
793: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
794: return(0);
795: }
799: /*@
800: TSThetaGetTheta - Get the abscissa of the stage in (0,1].
802: Not Collective
804: Input Parameter:
805: . ts - timestepping context
807: Output Parameter:
808: . theta - stage abscissa
810: Note:
811: Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.
813: Level: Advanced
815: .seealso: TSThetaSetTheta()
816: @*/
817: PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta)
818: {
824: PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
825: return(0);
826: }
830: /*@
831: TSThetaSetTheta - Set the abscissa of the stage in (0,1].
833: Not Collective
835: Input Parameter:
836: + ts - timestepping context
837: - theta - stage abscissa
839: Options Database:
840: . -ts_theta_theta <theta>
842: Level: Intermediate
844: .seealso: TSThetaGetTheta()
845: @*/
846: PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta)
847: {
852: PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
853: return(0);
854: }
858: /*@
859: TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
861: Not Collective
863: Input Parameter:
864: . ts - timestepping context
866: Output Parameter:
867: . endpoint - PETSC_TRUE when using the endpoint variant
869: Level: Advanced
871: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
872: @*/
873: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
874: {
880: PetscTryMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
881: return(0);
882: }
886: /*@
887: TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
889: Not Collective
891: Input Parameter:
892: + ts - timestepping context
893: - flg - PETSC_TRUE to use the endpoint variant
895: Options Database:
896: . -ts_theta_endpoint <flg>
898: Level: Intermediate
900: .seealso: TSTHETA, TSCN
901: @*/
902: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
903: {
908: PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
909: return(0);
910: }
912: /*
913: * TSBEULER and TSCN are straightforward specializations of TSTHETA.
914: * The creation functions for these specializations are below.
915: */
919: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
920: {
924: SNESView(ts->snes,viewer);
925: return(0);
926: }
928: /*MC
929: TSBEULER - ODE solver using the implicit backward Euler method
931: Level: beginner
933: Notes:
934: TSBEULER is equivalent to TSTHETA with Theta=1.0
936: $ -ts_type theta -ts_theta_theta 1.
938: .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA
940: M*/
943: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
944: {
948: TSCreate_Theta(ts);
949: TSThetaSetTheta(ts,1.0);
950: ts->ops->setup = TSSetUp_BEuler;
951: ts->ops->view = TSView_BEuler;
952: return(0);
953: }
957: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
958: {
962: SNESView(ts->snes,viewer);
963: return(0);
964: }
966: /*MC
967: TSCN - ODE solver using the implicit Crank-Nicolson method.
969: Level: beginner
971: Notes:
972: TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.
974: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint
976: .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA
978: M*/
981: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
982: {
986: TSCreate_Theta(ts);
987: TSThetaSetTheta(ts,0.5);
988: TSThetaSetEndpoint(ts,PETSC_TRUE);
989: ts->ops->setup = TSSetUp_CN;
990: ts->ops->view = TSView_CN;
991: return(0);
992: }