Actual source code: theta.c
petsc-3.8.4 2018-03-24
1: /*
2: Code for timestepping with implicit Theta method
3: */
4: #include <petsc/private/tsimpl.h>
5: #include <petscsnes.h>
6: #include <petscdm.h>
7: #include <petscmat.h>
9: typedef struct {
10: /* context for time stepping */
11: PetscReal stage_time;
12: Vec X0,X,Xdot; /* Storage for stages and time derivative */
13: Vec affine; /* Affine vector needed for residual at beginning of step in endpoint formulation */
14: PetscReal Theta;
15: PetscReal ptime;
16: PetscReal time_step;
17: PetscInt order;
18: PetscBool endpoint;
19: PetscBool extrapolate;
20: TSStepStatus status;
21: Vec VecCostIntegral0; /* Backup for roll-backs due to events */
23: /* context for sensitivity analysis */
24: PetscInt num_tlm; /* Total number of tangent linear equations */
25: Vec *VecsDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */
26: Vec *VecsDeltaMu; /* Increment of the adjoint sensitivity w.r.t P at stage */
27: Vec *VecsSensiTemp; /* Vector to be timed with Jacobian transpose */
28: Vec *VecsDeltaFwdSensi; /* Increment of the forward sensitivity at stage */
29: Vec VecIntegralSensiTemp; /* Working vector for forward integral sensitivity */
30: Vec VecIntegralSensipTemp; /* Working vector for forward integral sensitivity */
31: Vec *VecsFwdSensi0; /* backup for roll-backs due to events */
32: Vec *VecsIntegralSensi0; /* backup for roll-backs due to events */
33: Vec *VecsIntegralSensip0; /* backup for roll-backs due to events */
35: /* context for error estimation */
36: Vec vec_sol_prev;
37: Vec vec_lte_work;
38: } TS_Theta;
40: static PetscErrorCode TSThetaGetX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
41: {
42: TS_Theta *th = (TS_Theta*)ts->data;
46: if (X0) {
47: if (dm && dm != ts->dm) {
48: DMGetNamedGlobalVector(dm,"TSTheta_X0",X0);
49: } else *X0 = ts->vec_sol;
50: }
51: if (Xdot) {
52: if (dm && dm != ts->dm) {
53: DMGetNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
54: } else *Xdot = th->Xdot;
55: }
56: return(0);
57: }
59: static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
60: {
64: if (X0) {
65: if (dm && dm != ts->dm) {
66: DMRestoreNamedGlobalVector(dm,"TSTheta_X0",X0);
67: }
68: }
69: if (Xdot) {
70: if (dm && dm != ts->dm) {
71: DMRestoreNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
72: }
73: }
74: return(0);
75: }
77: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine,DM coarse,void *ctx)
78: {
80: return(0);
81: }
83: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
84: {
85: TS ts = (TS)ctx;
87: Vec X0,Xdot,X0_c,Xdot_c;
90: TSThetaGetX0AndXdot(ts,fine,&X0,&Xdot);
91: TSThetaGetX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
92: MatRestrict(restrct,X0,X0_c);
93: MatRestrict(restrct,Xdot,Xdot_c);
94: VecPointwiseMult(X0_c,rscale,X0_c);
95: VecPointwiseMult(Xdot_c,rscale,Xdot_c);
96: TSThetaRestoreX0AndXdot(ts,fine,&X0,&Xdot);
97: TSThetaRestoreX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
98: return(0);
99: }
101: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm,DM subdm,void *ctx)
102: {
104: return(0);
105: }
107: static PetscErrorCode DMSubDomainRestrictHook_TSTheta(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
108: {
109: TS ts = (TS)ctx;
111: Vec X0,Xdot,X0_sub,Xdot_sub;
114: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
115: TSThetaGetX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
117: VecScatterBegin(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
118: VecScatterEnd(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
120: VecScatterBegin(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
121: VecScatterEnd(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
123: TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
124: TSThetaRestoreX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
125: return(0);
126: }
128: static PetscErrorCode TSThetaEvaluateCostIntegral(TS ts)
129: {
130: TS_Theta *th = (TS_Theta*)ts->data;
134: if (th->endpoint) {
135: /* Evolve ts->vec_costintegral to compute integrals */
136: if (th->Theta!=1.0) {
137: TSComputeCostIntegrand(ts,th->ptime,th->X0,ts->vec_costintegrand);
138: VecAXPY(ts->vec_costintegral,th->time_step*(1.0-th->Theta),ts->vec_costintegrand);
139: }
140: TSComputeCostIntegrand(ts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
141: VecAXPY(ts->vec_costintegral,th->time_step*th->Theta,ts->vec_costintegrand);
142: } else {
143: TSComputeCostIntegrand(ts,th->stage_time,th->X,ts->vec_costintegrand);
144: VecAXPY(ts->vec_costintegral,th->time_step,ts->vec_costintegrand);
145: }
146: return(0);
147: }
149: static PetscErrorCode TSForwardCostIntegral_Theta(TS ts)
150: {
151: TS_Theta *th = (TS_Theta*)ts->data;
155: /* backup cost integral */
156: VecCopy(ts->vec_costintegral,th->VecCostIntegral0);
157: TSThetaEvaluateCostIntegral(ts);
158: return(0);
159: }
161: static PetscErrorCode TSAdjointCostIntegral_Theta(TS ts)
162: {
166: TSThetaEvaluateCostIntegral(ts);
167: return(0);
168: }
170: static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x)
171: {
172: PetscInt nits,lits;
176: SNESSolve(ts->snes,b,x);
177: SNESGetIterationNumber(ts->snes,&nits);
178: SNESGetLinearSolveIterations(ts->snes,&lits);
179: ts->snes_its += nits; ts->ksp_its += lits;
180: return(0);
181: }
183: static PetscErrorCode TSStep_Theta(TS ts)
184: {
185: TS_Theta *th = (TS_Theta*)ts->data;
186: PetscInt rejections = 0;
187: PetscBool stageok,accept = PETSC_TRUE;
188: PetscReal next_time_step = ts->time_step;
192: if (!ts->steprollback) {
193: if (th->vec_sol_prev) { VecCopy(th->X0,th->vec_sol_prev); }
194: VecCopy(ts->vec_sol,th->X0);
195: }
197: th->status = TS_STEP_INCOMPLETE;
198: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
200: PetscReal shift = 1/(th->Theta*ts->time_step);
201: th->stage_time = ts->ptime + (th->endpoint ? (PetscReal)1 : th->Theta)*ts->time_step;
203: VecCopy(th->X0,th->X);
204: if (th->extrapolate && !ts->steprestart) {
205: VecAXPY(th->X,1/shift,th->Xdot);
206: }
207: if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */
208: if (!th->affine) {VecDuplicate(ts->vec_sol,&th->affine);}
209: VecZeroEntries(th->Xdot);
210: TSComputeIFunction(ts,ts->ptime,th->X0,th->Xdot,th->affine,PETSC_FALSE);
211: VecScale(th->affine,(th->Theta-1)/th->Theta);
212: } else if (th->affine) { /* Just in case th->endpoint is changed between calls to TSStep_Theta() */
213: VecZeroEntries(th->affine);
214: }
215: TSPreStage(ts,th->stage_time);
216: TS_SNESSolve(ts,th->affine,th->X);
217: TSPostStage(ts,th->stage_time,0,&th->X);
218: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->X,&stageok);
219: if (!stageok) goto reject_step;
221: th->status = TS_STEP_PENDING;
222: if (th->endpoint) {
223: VecCopy(th->X,ts->vec_sol);
224: } else {
225: VecAXPBYPCZ(th->Xdot,-shift,shift,0,th->X0,th->X);
226: VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);
227: }
228: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
229: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
230: if (!accept) {
231: VecCopy(th->X0,ts->vec_sol);
232: ts->time_step = next_time_step;
233: goto reject_step;
234: }
236: if (ts->forward_solve || ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */
237: th->ptime = ts->ptime;
238: th->time_step = ts->time_step;
239: }
241: ts->ptime += ts->time_step;
242: ts->time_step = next_time_step;
243: break;
245: reject_step:
246: ts->reject++; accept = PETSC_FALSE;
247: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
248: ts->reason = TS_DIVERGED_STEP_REJECTED;
249: PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
250: }
251: }
252: return(0);
253: }
255: static PetscErrorCode TSAdjointStep_Theta(TS ts)
256: {
257: TS_Theta *th = (TS_Theta*)ts->data;
258: Vec *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
259: PetscInt nadj;
260: PetscErrorCode ierr;
261: Mat J,Jp;
262: KSP ksp;
263: PetscReal shift;
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 = th->endpoint ? ts->ptime : (ts->ptime+(1.-th->Theta)*ts->time_step); /* time_step is negative*/
272: th->ptime = ts->ptime + ts->time_step;
273: th->time_step = -ts->time_step;
275: /* Build RHS */
276: if (ts->vec_costintegral) { /* Cost function has an integral term */
277: if (th->endpoint) {
278: TSAdjointComputeDRDYFunction(ts,th->stage_time,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*th->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*th->time_step);
293: if (th->endpoint) {
294: TSComputeIJacobian(ts,th->stage_time,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.)*th->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: }
313: for (nadj=0; nadj<ts->numcost; nadj++) {
314: MatMultTranspose(J,VecsDeltaLam[nadj],ts->vecs_sensi[nadj]);
315: if (ts->vec_costintegral) {
316: VecAXPY(ts->vecs_sensi[nadj],-1.,ts->vecs_drdy[nadj]);
317: }
318: VecScale(ts->vecs_sensi[nadj],1./shift);
319: }
320: } else { /* backward Euler */
321: shift = 0.0;
322: TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
323: for (nadj=0; nadj<ts->numcost; nadj++) {
324: MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
325: VecAXPY(ts->vecs_sensi[nadj],-th->time_step,VecsSensiTemp[nadj]);
326: if (ts->vec_costintegral) {
327: VecAXPY(ts->vecs_sensi[nadj],th->time_step,ts->vecs_drdy[nadj]);
328: }
329: }
330: }
332: if (ts->vecs_sensip) { /* sensitivities wrt parameters */
333: TSAdjointComputeRHSJacobian(ts,th->stage_time,ts->vec_sol,ts->Jacp);
334: for (nadj=0; nadj<ts->numcost; nadj++) {
335: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
336: VecAXPY(ts->vecs_sensip[nadj],th->time_step*th->Theta,VecsDeltaMu[nadj]);
337: }
338: if (th->Theta!=1.) {
339: TSAdjointComputeRHSJacobian(ts,th->ptime,th->X0,ts->Jacp);
340: for (nadj=0; nadj<ts->numcost; nadj++) {
341: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
342: VecAXPY(ts->vecs_sensip[nadj],th->time_step*(1.-th->Theta),VecsDeltaMu[nadj]);
343: }
344: }
345: if (ts->vec_costintegral) {
346: TSAdjointComputeDRDPFunction(ts,th->stage_time,ts->vec_sol,ts->vecs_drdp);
347: for (nadj=0; nadj<ts->numcost; nadj++) {
348: VecAXPY(ts->vecs_sensip[nadj],th->time_step*th->Theta,ts->vecs_drdp[nadj]);
349: }
350: if (th->Theta!=1.) {
351: TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
352: for (nadj=0; nadj<ts->numcost; nadj++) {
353: VecAXPY(ts->vecs_sensip[nadj],th->time_step*(1.-th->Theta),ts->vecs_drdp[nadj]);
354: }
355: }
356: }
357: }
358: } else { /* one-stage case */
359: shift = 0.0;
360: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE); /* get -f_y */
361: if (ts->vec_costintegral) {
362: TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
363: }
364: for (nadj=0; nadj<ts->numcost; nadj++) {
365: MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
366: VecAXPY(ts->vecs_sensi[nadj],-th->time_step,VecsSensiTemp[nadj]);
367: if (ts->vec_costintegral) {
368: VecAXPY(ts->vecs_sensi[nadj],th->time_step,ts->vecs_drdy[nadj]);
369: }
370: }
371: if (ts->vecs_sensip) {
372: TSAdjointComputeRHSJacobian(ts,th->stage_time,th->X,ts->Jacp);
373: for (nadj=0; nadj<ts->numcost; nadj++) {
374: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
375: VecAXPY(ts->vecs_sensip[nadj],th->time_step,VecsDeltaMu[nadj]);
376: }
377: if (ts->vec_costintegral) {
378: TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
379: for (nadj=0; nadj<ts->numcost; nadj++) {
380: VecAXPY(ts->vecs_sensip[nadj],th->time_step,ts->vecs_drdp[nadj]);
381: }
382: }
383: }
384: }
386: th->status = TS_STEP_COMPLETE;
387: return(0);
388: }
390: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
391: {
392: TS_Theta *th = (TS_Theta*)ts->data;
393: PetscReal dt = t - ts->ptime;
397: VecCopy(ts->vec_sol,th->X);
398: if (th->endpoint) dt *= th->Theta;
399: VecWAXPY(X,dt,th->Xdot,th->X);
400: return(0);
401: }
403: static PetscErrorCode TSEvaluateWLTE_Theta(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
404: {
405: TS_Theta *th = (TS_Theta*)ts->data;
406: Vec X = ts->vec_sol; /* X = solution */
407: Vec Y = th->vec_lte_work; /* Y = X + LTE */
408: PetscReal wltea,wlter;
412: if (!th->vec_sol_prev) {*wlte = -1; return(0);}
413: /* Cannot compute LTE in first step or in restart after event */
414: if (ts->steprestart) {*wlte = -1; return(0);}
415: /* Compute LTE using backward differences with non-constant time step */
416: {
417: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
418: PetscReal a = 1 + h_prev/h;
419: PetscScalar scal[3]; Vec vecs[3];
420: scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
421: vecs[0] = X; vecs[1] = th->X0; vecs[2] = th->vec_sol_prev;
422: VecCopy(X,Y);
423: VecMAXPY(Y,3,scal,vecs);
424: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
425: }
426: if (order) *order = 2;
427: return(0);
428: }
430: static PetscErrorCode TSRollBack_Theta(TS ts)
431: {
432: TS_Theta *th = (TS_Theta*)ts->data;
433: PetscInt ntlm,ncost;
437: VecCopy(th->X0,ts->vec_sol);
438: if (ts->vec_costintegral && ts->costintegralfwd) {
439: VecCopy(th->VecCostIntegral0,ts->vec_costintegral);
440: }
441: th->status = TS_STEP_INCOMPLETE;
443: for (ntlm=0;ntlm<th->num_tlm;ntlm++) {
444: VecCopy(th->VecsFwdSensi0[ntlm],ts->vecs_fwdsensipacked[ntlm]);
445: }
446: if (ts->vecs_integral_sensi) {
447: for (ncost=0;ncost<ts->numcost;ncost++) {
448: VecCopy(th->VecsIntegralSensi0[ncost],ts->vecs_integral_sensi[ncost]);
449: }
450: }
451: if (ts->vecs_integral_sensip) {
452: for (ncost=0;ncost<ts->numcost;ncost++) {
453: VecCopy(th->VecsIntegralSensip0[ncost],ts->vecs_integral_sensip[ncost]);
454: }
455: }
456: return(0);
457: }
459: static PetscErrorCode TSThetaIntegrandDerivative(TS ts,PetscInt numtlm,Vec DRncostDY,Vec* s,Vec DRncostDP,Vec VecIntegrandDerivative)
460: {
461: PetscInt ntlm,low,high;
462: PetscScalar *v;
467: VecGetOwnershipRange(VecIntegrandDerivative,&low,&high);
468: VecGetArray(VecIntegrandDerivative,&v);
469: for (ntlm=low; ntlm<high; ntlm++) {
470: VecDot(DRncostDY,s[ntlm],&v[ntlm-low]);
471: }
472: VecRestoreArray(VecIntegrandDerivative,&v);
473: if (DRncostDP) {
474: VecAXPY(VecIntegrandDerivative,1.,DRncostDP);
475: }
476: return(0);
477: }
479: static PetscErrorCode TSForwardStep_Theta(TS ts)
480: {
481: TS_Theta *th = (TS_Theta*)ts->data;
482: Vec *VecsDeltaFwdSensi = th->VecsDeltaFwdSensi;
483: PetscInt ntlm,ncost,np;
484: KSP ksp;
485: Mat J,Jp;
486: PetscReal shift;
490: for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
491: VecCopy(ts->vecs_fwdsensipacked[ntlm],th->VecsFwdSensi0[ntlm]);
492: }
493: for (ncost=0; ncost<ts->numcost; ncost++) {
494: if (ts->vecs_integral_sensi) {
495: VecCopy(ts->vecs_integral_sensi[ncost],th->VecsIntegralSensi0[ncost]);
496: }
497: if (ts->vecs_integral_sensip) {
498: VecCopy(ts->vecs_integral_sensip[ncost],th->VecsIntegralSensip0[ncost]);
499: }
500: }
502: SNESGetKSP(ts->snes,&ksp);
503: TSGetIJacobian(ts,&J,&Jp,NULL,NULL);
505: /* Build RHS */
506: if (th->endpoint) { /* 2-stage method*/
507: shift = 1./((th->Theta-1.)*th->time_step);
508: TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,shift,J,Jp,PETSC_FALSE);
510: for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
511: MatMult(J,ts->vecs_fwdsensipacked[ntlm],VecsDeltaFwdSensi[ntlm]);
512: VecScale(VecsDeltaFwdSensi[ntlm],(th->Theta-1.)/th->Theta);
513: }
514: /* Add the f_p forcing terms */
515: TSForwardComputeRHSJacobianP(ts,th->ptime,th->X0,ts->vecs_jacp);
516: for (np=0; np<ts->num_parameters; np++) {
517: VecAXPY(VecsDeltaFwdSensi[np+ts->num_initialvalues],(1.-th->Theta)/th->Theta,ts->vecs_jacp[np]);
518: }
519: TSForwardComputeRHSJacobianP(ts,th->stage_time,ts->vec_sol,ts->vecs_jacp);
520: for (np=0; np<ts->num_parameters; np++) {
521: VecAXPY(VecsDeltaFwdSensi[np+ts->num_initialvalues],1,ts->vecs_jacp[np]);
522: }
523: } else { /* 1-stage method */
524: shift = 0.0;
525: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
526: for (ntlm=0; ntlm<ts->num_parameters; ntlm++) {
527: MatMult(J,ts->vecs_fwdsensipacked[ntlm],VecsDeltaFwdSensi[ntlm]);
528: VecScale(VecsDeltaFwdSensi[ntlm],-1);
529: }
530: /* Add the f_p forcing terms */
531: TSForwardComputeRHSJacobianP(ts,th->stage_time,th->X,ts->vecs_jacp);
532: for (np=0; np<ts->num_parameters; np++) {
533: VecAXPY(VecsDeltaFwdSensi[np+ts->num_initialvalues],1,ts->vecs_jacp[np]);
534: }
535: }
537: /* Build LHS */
538: shift = 1/(th->Theta*th->time_step);
539: if (th->endpoint) {
540: TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE);
541: } else {
542: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
543: }
544: KSPSetOperators(ksp,J,Jp);
546: /*
547: Evaluate the first stage of integral gradients with the 2-stage method:
548: drdy|t_n*S(t_n) + drdp|t_n
549: This is done before the linear solve because the sensitivity variable S(t_n) will be propagated to S(t_{n+1})
550: It is important to note that sensitivitesi to parameters (sensip) and sensitivities to initial values (sensi) are independent of each other, but integral sensip relies on sensip and integral sensi requires sensi
551: */
552: if (th->endpoint) { /* 2-stage method only */
553: if (ts->vecs_integral_sensi || ts->vecs_integral_sensip) {
554: TSAdjointComputeDRDYFunction(ts,th->ptime,th->X0,ts->vecs_drdy);
555: }
556: if (ts->vecs_integral_sensip) {
557: TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
558: for (ncost=0; ncost<ts->numcost; ncost++) {
559: TSThetaIntegrandDerivative(ts,ts->num_parameters,ts->vecs_drdy[ncost],&ts->vecs_fwdsensipacked[ts->num_initialvalues],ts->vecs_drdp[ncost],th->VecIntegralSensipTemp);
560: VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*(1.-th->Theta),th->VecIntegralSensipTemp);
561: }
562: }
563: if (ts->vecs_integral_sensi) {
564: for (ncost=0; ncost<ts->numcost; ncost++) {
565: TSThetaIntegrandDerivative(ts,ts->num_initialvalues,ts->vecs_drdy[ncost],ts->vecs_fwdsensipacked,NULL,th->VecIntegralSensiTemp);
566: VecAXPY(ts->vecs_integral_sensi[ncost],th->time_step*(1.-th->Theta),th->VecIntegralSensiTemp);
567: }
568: }
569: }
571: /* Solve the tangent linear equation for forward sensitivities to parameters */
572: for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
573: if (th->endpoint) {
574: KSPSolve(ksp,VecsDeltaFwdSensi[ntlm],ts->vecs_fwdsensipacked[ntlm]);
575: } else {
576: KSPSolve(ksp,VecsDeltaFwdSensi[ntlm],VecsDeltaFwdSensi[ntlm]);
577: VecAXPY(ts->vecs_fwdsensipacked[ntlm],1.,VecsDeltaFwdSensi[ntlm]);
578: }
579: }
580: /* Evaluate the second stage of integral gradients with the 2-stage method:
581: drdy|t_{n+1}*S(t_{n+1}) + drdp|t_{n+1}
582: */
583: if (ts->vecs_integral_sensi || ts->vecs_integral_sensip) {
584: TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
585: }
586: if (ts->vecs_integral_sensip) {
587: TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
588: for (ncost=0; ncost<ts->numcost; ncost++) {
589: TSThetaIntegrandDerivative(ts,ts->num_parameters,ts->vecs_drdy[ncost],&ts->vecs_fwdsensipacked[ts->num_initialvalues],ts->vecs_drdp[ncost],th->VecIntegralSensipTemp);
590: if (th->endpoint) {
591: VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*th->Theta,th->VecIntegralSensipTemp);
592: } else {
593: VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step,th->VecIntegralSensipTemp);
594: }
595: }
596: }
597: if (ts->vecs_integral_sensi) {
598: for (ncost=0; ncost<ts->numcost; ncost++) {
599: TSThetaIntegrandDerivative(ts,ts->num_initialvalues,ts->vecs_drdy[ncost],ts->vecs_fwdsensipacked,NULL,th->VecIntegralSensiTemp);
600: if (th->endpoint) {
601: VecAXPY(ts->vecs_integral_sensi[ncost],th->time_step*th->Theta,th->VecIntegralSensiTemp);
602: } else {
603: VecAXPY(ts->vecs_integral_sensi[ncost],th->time_step,th->VecIntegralSensiTemp);
604: }
605: }
606: }
608: if (!th->endpoint) { /* VecsDeltaFwdSensip are the increment for the stage values for the 1-stage method */
609: for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
610: VecAXPY(ts->vecs_fwdsensipacked[ntlm],(1.-th->Theta)/th->Theta,VecsDeltaFwdSensi[ntlm]);
611: }
612: }
613: return(0);
614: }
616: /*------------------------------------------------------------*/
617: static PetscErrorCode TSReset_Theta(TS ts)
618: {
619: TS_Theta *th = (TS_Theta*)ts->data;
623: VecDestroy(&th->X);
624: VecDestroy(&th->Xdot);
625: VecDestroy(&th->X0);
626: VecDestroy(&th->affine);
628: VecDestroy(&th->vec_sol_prev);
629: VecDestroy(&th->vec_lte_work);
631: VecDestroy(&th->VecCostIntegral0);
632: if (ts->forward_solve) {
633: VecDestroyVecs(th->num_tlm,&th->VecsDeltaFwdSensi);
634: VecDestroyVecs(th->num_tlm,&th->VecsFwdSensi0);
635: if (ts->vecs_integral_sensi) {
636: VecDestroy(&th->VecIntegralSensiTemp);
637: VecDestroyVecs(ts->numcost,&th->VecsIntegralSensi0);
638: }
639: if (ts->vecs_integral_sensip) {
640: VecDestroy(&th->VecIntegralSensipTemp);
641: VecDestroyVecs(ts->numcost,&th->VecsIntegralSensip0);
642: }
643: }
644: VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
645: VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
646: VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);
648: return(0);
649: }
651: static PetscErrorCode TSDestroy_Theta(TS ts)
652: {
656: TSReset_Theta(ts);
657: if (ts->dm) {
658: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
659: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
660: }
661: PetscFree(ts->data);
662: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
663: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
664: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
665: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
666: return(0);
667: }
669: /*
670: This defines the nonlinear equation that is to be solved with SNES
671: G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
672: */
673: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
674: {
675: TS_Theta *th = (TS_Theta*)ts->data;
677: Vec X0,Xdot;
678: DM dm,dmsave;
679: PetscReal shift = 1/(th->Theta*ts->time_step);
682: SNESGetDM(snes,&dm);
683: /* When using the endpoint variant, this is actually 1/Theta * Xdot */
684: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
685: VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);
687: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
688: dmsave = ts->dm;
689: ts->dm = dm;
690: TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
691: ts->dm = dmsave;
692: TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
693: return(0);
694: }
696: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
697: {
698: TS_Theta *th = (TS_Theta*)ts->data;
700: Vec Xdot;
701: DM dm,dmsave;
702: PetscReal shift = 1/(th->Theta*ts->time_step);
705: SNESGetDM(snes,&dm);
706: /* Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
707: TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);
709: dmsave = ts->dm;
710: ts->dm = dm;
711: TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
712: ts->dm = dmsave;
713: TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
714: return(0);
715: }
717: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
718: {
719: TS_Theta *th = (TS_Theta*)ts->data;
723: /* combine sensitivities to parameters and sensitivities to initial values into one array */
724: th->num_tlm = ts->num_parameters + ts->num_initialvalues;
725: VecDuplicateVecs(ts->vecs_fwdsensipacked[0],th->num_tlm,&th->VecsDeltaFwdSensi);
726: if (ts->vecs_integral_sensi) {
727: VecDuplicate(ts->vecs_integral_sensi[0],&th->VecIntegralSensiTemp);
728: }
729: if (ts->vecs_integral_sensip) {
730: VecDuplicate(ts->vecs_integral_sensip[0],&th->VecIntegralSensipTemp);
731: }
732: /* backup sensitivity results for roll-backs */
733: VecDuplicateVecs(ts->vecs_fwdsensipacked[0],th->num_tlm,&th->VecsFwdSensi0);
734: if (ts->vecs_integral_sensi) {
735: VecDuplicateVecs(ts->vecs_integral_sensi[0],ts->numcost,&th->VecsIntegralSensi0);
736: }
737: if (ts->vecs_integral_sensip) {
738: VecDuplicateVecs(ts->vecs_integral_sensip[0],ts->numcost,&th->VecsIntegralSensip0);
739: }
740: return(0);
741: }
743: static PetscErrorCode TSSetUp_Theta(TS ts)
744: {
745: TS_Theta *th = (TS_Theta*)ts->data;
746: PetscBool match;
750: if (!th->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */
751: VecDuplicate(ts->vec_costintegral,&th->VecCostIntegral0);
752: }
753: if (!th->X) {
754: VecDuplicate(ts->vec_sol,&th->X);
755: }
756: if (!th->Xdot) {
757: VecDuplicate(ts->vec_sol,&th->Xdot);
758: }
759: if (!th->X0) {
760: VecDuplicate(ts->vec_sol,&th->X0);
761: }
762: if (th->endpoint) {
763: VecDuplicate(ts->vec_sol,&th->affine);
764: }
766: th->order = (th->Theta == 0.5) ? 2 : 1;
768: TSGetDM(ts,&ts->dm);
769: DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
770: DMSubDomainHookAdd(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
772: TSGetAdapt(ts,&ts->adapt);
773: TSAdaptCandidatesClear(ts->adapt);
774: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
775: if (!match) {
776: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
777: VecDuplicate(ts->vec_sol,&th->vec_lte_work);
778: }
779: TSGetSNES(ts,&ts->snes);
780: return(0);
781: }
783: /*------------------------------------------------------------*/
785: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
786: {
787: TS_Theta *th = (TS_Theta*)ts->data;
791: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
792: if(ts->vecs_sensip) {
793: VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
794: }
795: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
796: return(0);
797: }
799: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
800: {
801: TS_Theta *th = (TS_Theta*)ts->data;
805: PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
806: {
807: PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
808: PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
809: PetscOptionsBool("-ts_theta_initial_guess_extrapolate","Extrapolate stage initial guess from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
810: }
811: PetscOptionsTail();
812: return(0);
813: }
815: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
816: {
817: TS_Theta *th = (TS_Theta*)ts->data;
818: PetscBool iascii;
822: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
823: if (iascii) {
824: PetscViewerASCIIPrintf(viewer," Theta=%g\n",(double)th->Theta);
825: PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
826: }
827: return(0);
828: }
830: static PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
831: {
832: TS_Theta *th = (TS_Theta*)ts->data;
835: *theta = th->Theta;
836: return(0);
837: }
839: static PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
840: {
841: TS_Theta *th = (TS_Theta*)ts->data;
844: if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
845: th->Theta = theta;
846: th->order = (th->Theta == 0.5) ? 2 : 1;
847: return(0);
848: }
850: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
851: {
852: TS_Theta *th = (TS_Theta*)ts->data;
855: *endpoint = th->endpoint;
856: return(0);
857: }
859: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
860: {
861: TS_Theta *th = (TS_Theta*)ts->data;
864: th->endpoint = flg;
865: return(0);
866: }
868: #if defined(PETSC_HAVE_COMPLEX)
869: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
870: {
871: PetscComplex z = xr + xi*PETSC_i,f;
872: TS_Theta *th = (TS_Theta*)ts->data;
873: const PetscReal one = 1.0;
876: f = (one + (one - th->Theta)*z)/(one - th->Theta*z);
877: *yr = PetscRealPartComplex(f);
878: *yi = PetscImaginaryPartComplex(f);
879: return(0);
880: }
881: #endif
883: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
884: {
885: TS_Theta *th = (TS_Theta*)ts->data;
888: if (ns) *ns = 1;
889: if (Y) *Y = th->endpoint ? &(th->X0) : &(th->X);
890: return(0);
891: }
893: /* ------------------------------------------------------------ */
894: /*MC
895: TSTHETA - DAE solver using the implicit Theta method
897: Level: beginner
899: Options Database:
900: + -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
901: . -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of the Theta method
902: - -ts_theta_initial_guess_extrapolate <flg> - Extrapolate stage initial guess from previous solution (sometimes unstable)
904: Notes:
905: $ -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
906: $ -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
907: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)
909: This method can be applied to DAE.
911: This method is cast as a 1-stage implicit Runge-Kutta method.
913: .vb
914: Theta | Theta
915: -------------
916: | 1
917: .ve
919: For the default Theta=0.5, this is also known as the implicit midpoint rule.
921: When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:
923: .vb
924: 0 | 0 0
925: 1 | 1-Theta Theta
926: -------------------
927: | 1-Theta Theta
928: .ve
930: For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).
932: To apply a diagonally implicit RK method to DAE, the stage formula
934: $ Y_i = X + h sum_j a_ij Y'_j
936: is interpreted as a formula for Y'_i in terms of Y_i and known values (Y'_j, j<i)
938: .seealso: TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()
940: M*/
941: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
942: {
943: TS_Theta *th;
947: ts->ops->reset = TSReset_Theta;
948: ts->ops->destroy = TSDestroy_Theta;
949: ts->ops->view = TSView_Theta;
950: ts->ops->setup = TSSetUp_Theta;
951: ts->ops->adjointsetup = TSAdjointSetUp_Theta;
952: ts->ops->step = TSStep_Theta;
953: ts->ops->interpolate = TSInterpolate_Theta;
954: ts->ops->evaluatewlte = TSEvaluateWLTE_Theta;
955: ts->ops->rollback = TSRollBack_Theta;
956: ts->ops->setfromoptions = TSSetFromOptions_Theta;
957: ts->ops->snesfunction = SNESTSFormFunction_Theta;
958: ts->ops->snesjacobian = SNESTSFormJacobian_Theta;
959: #if defined(PETSC_HAVE_COMPLEX)
960: ts->ops->linearstability = TSComputeLinearStability_Theta;
961: #endif
962: ts->ops->getstages = TSGetStages_Theta;
963: ts->ops->adjointstep = TSAdjointStep_Theta;
964: ts->ops->adjointintegral = TSAdjointCostIntegral_Theta;
965: ts->ops->forwardintegral = TSForwardCostIntegral_Theta;
966: ts->default_adapt_type = TSADAPTNONE;
967: ts->ops->forwardsetup = TSForwardSetUp_Theta;
968: ts->ops->forwardstep = TSForwardStep_Theta;
970: ts->usessnes = PETSC_TRUE;
972: PetscNewLog(ts,&th);
973: ts->data = (void*)th;
975: th->VecsDeltaLam = NULL;
976: th->VecsDeltaMu = NULL;
977: th->VecsSensiTemp = NULL;
979: th->extrapolate = PETSC_FALSE;
980: th->Theta = 0.5;
981: th->order = 2;
982: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
983: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
984: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
985: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
986: return(0);
987: }
989: /*@
990: TSThetaGetTheta - Get the abscissa of the stage in (0,1].
992: Not Collective
994: Input Parameter:
995: . ts - timestepping context
997: Output Parameter:
998: . theta - stage abscissa
1000: Note:
1001: Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.
1003: Level: Advanced
1005: .seealso: TSThetaSetTheta()
1006: @*/
1007: PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta)
1008: {
1014: PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
1015: return(0);
1016: }
1018: /*@
1019: TSThetaSetTheta - Set the abscissa of the stage in (0,1].
1021: Not Collective
1023: Input Parameter:
1024: + ts - timestepping context
1025: - theta - stage abscissa
1027: Options Database:
1028: . -ts_theta_theta <theta>
1030: Level: Intermediate
1032: .seealso: TSThetaGetTheta()
1033: @*/
1034: PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta)
1035: {
1040: PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
1041: return(0);
1042: }
1044: /*@
1045: TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
1047: Not Collective
1049: Input Parameter:
1050: . ts - timestepping context
1052: Output Parameter:
1053: . endpoint - PETSC_TRUE when using the endpoint variant
1055: Level: Advanced
1057: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
1058: @*/
1059: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
1060: {
1066: PetscUseMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
1067: return(0);
1068: }
1070: /*@
1071: TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
1073: Not Collective
1075: Input Parameter:
1076: + ts - timestepping context
1077: - flg - PETSC_TRUE to use the endpoint variant
1079: Options Database:
1080: . -ts_theta_endpoint <flg>
1082: Level: Intermediate
1084: .seealso: TSTHETA, TSCN
1085: @*/
1086: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
1087: {
1092: PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
1093: return(0);
1094: }
1096: /*
1097: * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1098: * The creation functions for these specializations are below.
1099: */
1101: static PetscErrorCode TSSetUp_BEuler(TS ts)
1102: {
1103: TS_Theta *th = (TS_Theta*)ts->data;
1107: 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");
1108: if (th->endpoint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_OPT_OVERWRITE,"Can not change to the endpoint form of the Theta methods when using backward Euler\n");
1109: TSSetUp_Theta(ts);
1110: return(0);
1111: }
1113: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
1114: {
1116: return(0);
1117: }
1119: /*MC
1120: TSBEULER - ODE solver using the implicit backward Euler method
1122: Level: beginner
1124: Notes:
1125: TSBEULER is equivalent to TSTHETA with Theta=1.0
1127: $ -ts_type theta -ts_theta_theta 1.0
1129: .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA
1131: M*/
1132: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1133: {
1137: TSCreate_Theta(ts);
1138: TSThetaSetTheta(ts,1.0);
1139: TSThetaSetEndpoint(ts,PETSC_FALSE);
1140: ts->ops->setup = TSSetUp_BEuler;
1141: ts->ops->view = TSView_BEuler;
1142: return(0);
1143: }
1145: static PetscErrorCode TSSetUp_CN(TS ts)
1146: {
1147: TS_Theta *th = (TS_Theta*)ts->data;
1151: 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");
1152: 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");
1153: TSSetUp_Theta(ts);
1154: return(0);
1155: }
1157: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
1158: {
1160: return(0);
1161: }
1163: /*MC
1164: TSCN - ODE solver using the implicit Crank-Nicolson method.
1166: Level: beginner
1168: Notes:
1169: TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.
1171: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint
1173: .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA
1175: M*/
1176: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1177: {
1181: TSCreate_Theta(ts);
1182: TSThetaSetTheta(ts,0.5);
1183: TSThetaSetEndpoint(ts,PETSC_TRUE);
1184: ts->ops->setup = TSSetUp_CN;
1185: ts->ops->view = TSView_CN;
1186: return(0);
1187: }