Actual source code: theta.c
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 Stages[2]; /* Storage for stage solutions */
13: Vec X0,X,Xdot; /* Storage for u^n, u^n + dt a_{11} k_1, and time derivative u^{n+1}_t */
14: Vec affine; /* Affine vector needed for residual at beginning of step in endpoint formulation */
15: PetscReal Theta;
16: PetscReal shift; /* Shift parameter for SNES Jacobian, used by forward, TLM and adjoint */
17: PetscInt order;
18: PetscBool endpoint;
19: PetscBool extrapolate;
20: TSStepStatus status;
21: Vec VecCostIntegral0; /* Backup for roll-backs due to events, used by cost integral */
22: PetscReal ptime0; /* Backup for ts->ptime, the start time of current time step, used by TLM and cost integral */
23: PetscReal time_step0; /* Backup for ts->timestep, the step size of current time step, used by TLM and cost integral*/
25: /* context for sensitivity analysis */
26: PetscInt num_tlm; /* Total number of tangent linear equations */
27: Vec *VecsDeltaLam; /* Increment of the adjoint sensitivity w.r.t IC at stage */
28: Vec *VecsDeltaMu; /* Increment of the adjoint sensitivity w.r.t P at stage */
29: Vec *VecsSensiTemp; /* Vector to be multiplied with Jacobian transpose */
30: Mat MatFwdStages[2]; /* TLM Stages */
31: Mat MatDeltaFwdSensip; /* Increment of the forward sensitivity at stage */
32: Vec VecDeltaFwdSensipCol; /* Working vector for holding one column of the sensitivity matrix */
33: Mat MatFwdSensip0; /* backup for roll-backs due to events */
34: Mat MatIntegralSensipTemp; /* Working vector for forward integral sensitivity */
35: Mat MatIntegralSensip0; /* backup for roll-backs due to events */
36: Vec *VecsDeltaLam2; /* Increment of the 2nd-order adjoint sensitivity w.r.t IC at stage */
37: Vec *VecsDeltaMu2; /* Increment of the 2nd-order adjoint sensitivity w.r.t P at stage */
38: Vec *VecsSensi2Temp; /* Working vectors that holds the residual for the second-order adjoint */
39: Vec *VecsAffine; /* Working vectors to store residuals */
40: /* context for error estimation */
41: Vec vec_sol_prev;
42: Vec vec_lte_work;
43: } TS_Theta;
45: static PetscErrorCode TSThetaGetX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
46: {
47: TS_Theta *th = (TS_Theta*)ts->data;
49: if (X0) {
50: if (dm && dm != ts->dm) {
51: DMGetNamedGlobalVector(dm,"TSTheta_X0",X0);
52: } else *X0 = ts->vec_sol;
53: }
54: if (Xdot) {
55: if (dm && dm != ts->dm) {
56: DMGetNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
57: } else *Xdot = th->Xdot;
58: }
59: return 0;
60: }
62: static PetscErrorCode TSThetaRestoreX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
63: {
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: {
79: return 0;
80: }
82: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
83: {
84: TS ts = (TS)ctx;
85: 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: }
98: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm,DM subdm,void *ctx)
99: {
100: return 0;
101: }
103: static PetscErrorCode DMSubDomainRestrictHook_TSTheta(DM dm,VecScatter gscat,VecScatter lscat,DM subdm,void *ctx)
104: {
105: TS ts = (TS)ctx;
106: Vec X0,Xdot,X0_sub,Xdot_sub;
108: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
109: TSThetaGetX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
111: VecScatterBegin(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
112: VecScatterEnd(gscat,X0,X0_sub,INSERT_VALUES,SCATTER_FORWARD);
114: VecScatterBegin(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
115: VecScatterEnd(gscat,Xdot,Xdot_sub,INSERT_VALUES,SCATTER_FORWARD);
117: TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
118: TSThetaRestoreX0AndXdot(ts,subdm,&X0_sub,&Xdot_sub);
119: return 0;
120: }
122: static PetscErrorCode TSThetaEvaluateCostIntegral(TS ts)
123: {
124: TS_Theta *th = (TS_Theta*)ts->data;
125: TS quadts = ts->quadraturets;
127: if (th->endpoint) {
128: /* Evolve ts->vec_costintegral to compute integrals */
129: if (th->Theta!=1.0) {
130: TSComputeRHSFunction(quadts,th->ptime0,th->X0,ts->vec_costintegrand);
131: VecAXPY(quadts->vec_sol,th->time_step0*(1.0-th->Theta),ts->vec_costintegrand);
132: }
133: TSComputeRHSFunction(quadts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
134: VecAXPY(quadts->vec_sol,th->time_step0*th->Theta,ts->vec_costintegrand);
135: } else {
136: TSComputeRHSFunction(quadts,th->stage_time,th->X,ts->vec_costintegrand);
137: VecAXPY(quadts->vec_sol,th->time_step0,ts->vec_costintegrand);
138: }
139: return 0;
140: }
142: static PetscErrorCode TSForwardCostIntegral_Theta(TS ts)
143: {
144: TS_Theta *th = (TS_Theta*)ts->data;
145: TS quadts = ts->quadraturets;
147: /* backup cost integral */
148: VecCopy(quadts->vec_sol,th->VecCostIntegral0);
149: TSThetaEvaluateCostIntegral(ts);
150: return 0;
151: }
153: static PetscErrorCode TSAdjointCostIntegral_Theta(TS ts)
154: {
155: TS_Theta *th = (TS_Theta*)ts->data;
157: /* Like TSForwardCostIntegral(), the adjoint cost integral evaluation relies on ptime0 and time_step0. */
158: th->ptime0 = ts->ptime + ts->time_step;
159: th->time_step0 = -ts->time_step;
160: TSThetaEvaluateCostIntegral(ts);
161: return 0;
162: }
164: static PetscErrorCode TSTheta_SNESSolve(TS ts,Vec b,Vec x)
165: {
166: PetscInt nits,lits;
168: SNESSolve(ts->snes,b,x);
169: SNESGetIterationNumber(ts->snes,&nits);
170: SNESGetLinearSolveIterations(ts->snes,&lits);
171: ts->snes_its += nits; ts->ksp_its += lits;
172: return 0;
173: }
175: static PetscErrorCode TSStep_Theta(TS ts)
176: {
177: TS_Theta *th = (TS_Theta*)ts->data;
178: PetscInt rejections = 0;
179: PetscBool stageok,accept = PETSC_TRUE;
180: PetscReal next_time_step = ts->time_step;
182: if (!ts->steprollback) {
183: if (th->vec_sol_prev) VecCopy(th->X0,th->vec_sol_prev);
184: VecCopy(ts->vec_sol,th->X0);
185: }
187: th->status = TS_STEP_INCOMPLETE;
188: while (!ts->reason && th->status != TS_STEP_COMPLETE) {
189: th->shift = 1/(th->Theta*ts->time_step);
190: th->stage_time = ts->ptime + (th->endpoint ? (PetscReal)1 : th->Theta)*ts->time_step;
191: VecCopy(th->X0,th->X);
192: if (th->extrapolate && !ts->steprestart) {
193: VecAXPY(th->X,1/th->shift,th->Xdot);
194: }
195: if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */
196: if (!th->affine) VecDuplicate(ts->vec_sol,&th->affine);
197: VecZeroEntries(th->Xdot);
198: TSComputeIFunction(ts,ts->ptime,th->X0,th->Xdot,th->affine,PETSC_FALSE);
199: VecScale(th->affine,(th->Theta-1)/th->Theta);
200: } else if (th->affine) { /* Just in case th->endpoint is changed between calls to TSStep_Theta() */
201: VecZeroEntries(th->affine);
202: }
203: TSPreStage(ts,th->stage_time);
204: TSTheta_SNESSolve(ts,th->affine,th->X);
205: TSPostStage(ts,th->stage_time,0,&th->X);
206: TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->X,&stageok);
207: if (!stageok) goto reject_step;
209: th->status = TS_STEP_PENDING;
210: if (th->endpoint) {
211: VecCopy(th->X,ts->vec_sol);
212: } else {
213: VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->X0,th->X);
214: VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);
215: }
216: TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
217: th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
218: if (!accept) {
219: VecCopy(th->X0,ts->vec_sol);
220: ts->time_step = next_time_step;
221: goto reject_step;
222: }
224: if (ts->forward_solve || ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */
225: th->ptime0 = ts->ptime;
226: th->time_step0 = ts->time_step;
227: }
228: ts->ptime += ts->time_step;
229: ts->time_step = next_time_step;
230: break;
232: reject_step:
233: ts->reject++; accept = PETSC_FALSE;
234: if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
235: ts->reason = TS_DIVERGED_STEP_REJECTED;
236: PetscInfo(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
237: }
238: }
239: return 0;
240: }
242: static PetscErrorCode TSAdjointStepBEuler_Private(TS ts)
243: {
244: TS_Theta *th = (TS_Theta*)ts->data;
245: TS quadts = ts->quadraturets;
246: Vec *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
247: Vec *VecsDeltaLam2 = th->VecsDeltaLam2,*VecsDeltaMu2 = th->VecsDeltaMu2,*VecsSensi2Temp = th->VecsSensi2Temp;
248: PetscInt nadj;
249: Mat J,Jpre,quadJ = NULL,quadJp = NULL;
250: KSP ksp;
251: PetscScalar *xarr;
252: TSEquationType eqtype;
253: PetscBool isexplicitode = PETSC_FALSE;
254: PetscReal adjoint_time_step;
256: TSGetEquationType(ts,&eqtype);
257: if (eqtype == TS_EQ_ODE_EXPLICIT) {
258: isexplicitode = PETSC_TRUE;
259: VecsDeltaLam = ts->vecs_sensi;
260: VecsDeltaLam2 = ts->vecs_sensi2;
261: }
262: th->status = TS_STEP_INCOMPLETE;
263: SNESGetKSP(ts->snes,&ksp);
264: TSGetIJacobian(ts,&J,&Jpre,NULL,NULL);
265: if (quadts) {
266: TSGetRHSJacobian(quadts,&quadJ,NULL,NULL,NULL);
267: TSGetRHSJacobianP(quadts,&quadJp,NULL,NULL);
268: }
270: th->stage_time = ts->ptime;
271: adjoint_time_step = -ts->time_step; /* always positive since time_step is negative */
273: /* Build RHS for first-order adjoint lambda_{n+1}/h + r_u^T(n+1) */
274: if (quadts) {
275: TSComputeRHSJacobian(quadts,th->stage_time,ts->vec_sol,quadJ,NULL);
276: }
278: for (nadj=0; nadj<ts->numcost; nadj++) {
279: VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);
280: VecScale(VecsSensiTemp[nadj],1./adjoint_time_step); /* lambda_{n+1}/h */
281: if (quadJ) {
282: MatDenseGetColumn(quadJ,nadj,&xarr);
283: VecPlaceArray(ts->vec_drdu_col,xarr);
284: VecAXPY(VecsSensiTemp[nadj],1.,ts->vec_drdu_col);
285: VecResetArray(ts->vec_drdu_col);
286: MatDenseRestoreColumn(quadJ,&xarr);
287: }
288: }
290: /* Build LHS for first-order adjoint */
291: th->shift = 1./adjoint_time_step;
292: TSComputeSNESJacobian(ts,ts->vec_sol,J,Jpre);
293: KSPSetOperators(ksp,J,Jpre);
295: /* Solve stage equation LHS*lambda_s = RHS for first-order adjoint */
296: for (nadj=0; nadj<ts->numcost; nadj++) {
297: KSPConvergedReason kspreason;
298: KSPSolveTranspose(ksp,VecsSensiTemp[nadj],VecsDeltaLam[nadj]);
299: KSPGetConvergedReason(ksp,&kspreason);
300: if (kspreason < 0) {
301: ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
302: PetscInfo(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n",ts->steps,nadj);
303: }
304: }
306: if (ts->vecs_sensi2) { /* U_{n+1} */
307: /* Get w1 at t_{n+1} from TLM matrix */
308: MatDenseGetColumn(ts->mat_sensip,0,&xarr);
309: VecPlaceArray(ts->vec_sensip_col,xarr);
310: /* lambda_s^T F_UU w_1 */
311: TSComputeIHessianProductFunctionUU(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fuu);
312: /* lambda_s^T F_UP w_2 */
313: TSComputeIHessianProductFunctionUP(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_dir,ts->vecs_fup);
314: for (nadj=0; nadj<ts->numcost; nadj++) { /* compute the residual */
315: VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);
316: VecScale(VecsSensi2Temp[nadj],1./adjoint_time_step);
317: VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fuu[nadj]);
318: if (ts->vecs_fup) {
319: VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fup[nadj]);
320: }
321: }
322: /* Solve stage equation LHS X = RHS for second-order adjoint */
323: for (nadj=0; nadj<ts->numcost; nadj++) {
324: KSPConvergedReason kspreason;
325: KSPSolveTranspose(ksp,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj]);
326: KSPGetConvergedReason(ksp,&kspreason);
327: if (kspreason < 0) {
328: ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
329: PetscInfo(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 2nd-order adjoint solve\n",ts->steps,nadj);
330: }
331: }
332: }
334: /* Update sensitivities, and evaluate integrals if there is any */
335: if (!isexplicitode) {
336: th->shift = 0.0;
337: TSComputeSNESJacobian(ts,ts->vec_sol,J,Jpre);
338: KSPSetOperators(ksp,J,Jpre);
339: for (nadj=0; nadj<ts->numcost; nadj++) {
340: /* Add f_U \lambda_s to the original RHS */
341: VecScale(VecsSensiTemp[nadj],-1.);
342: MatMultTransposeAdd(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj],VecsSensiTemp[nadj]);
343: VecScale(VecsSensiTemp[nadj],-adjoint_time_step);
344: VecCopy(VecsSensiTemp[nadj],ts->vecs_sensi[nadj]);
345: if (ts->vecs_sensi2) {
346: MatMultTransposeAdd(J,VecsDeltaLam2[nadj],VecsSensi2Temp[nadj],VecsSensi2Temp[nadj]);
347: VecScale(VecsSensi2Temp[nadj],-adjoint_time_step);
348: VecCopy(VecsSensi2Temp[nadj],ts->vecs_sensi2[nadj]);
349: }
350: }
351: }
352: if (ts->vecs_sensip) {
353: VecAXPBYPCZ(th->Xdot,-1./adjoint_time_step,1.0/adjoint_time_step,0,th->X0,ts->vec_sol);
354: TSComputeIJacobianP(ts,th->stage_time,ts->vec_sol,th->Xdot,1./adjoint_time_step,ts->Jacp,PETSC_FALSE); /* get -f_p */
355: if (quadts) {
356: TSComputeRHSJacobianP(quadts,th->stage_time,ts->vec_sol,quadJp);
357: }
358: if (ts->vecs_sensi2p) {
359: /* lambda_s^T F_PU w_1 */
360: TSComputeIHessianProductFunctionPU(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fpu);
361: /* lambda_s^T F_PP w_2 */
362: TSComputeIHessianProductFunctionPP(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_dir,ts->vecs_fpp);
363: }
365: for (nadj=0; nadj<ts->numcost; nadj++) {
366: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
367: VecAXPY(ts->vecs_sensip[nadj],-adjoint_time_step,VecsDeltaMu[nadj]);
368: if (quadJp) {
369: MatDenseGetColumn(quadJp,nadj,&xarr);
370: VecPlaceArray(ts->vec_drdp_col,xarr);
371: VecAXPY(ts->vecs_sensip[nadj],adjoint_time_step,ts->vec_drdp_col);
372: VecResetArray(ts->vec_drdp_col);
373: MatDenseRestoreColumn(quadJp,&xarr);
374: }
375: if (ts->vecs_sensi2p) {
376: MatMultTranspose(ts->Jacp,VecsDeltaLam2[nadj],VecsDeltaMu2[nadj]);
377: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step,VecsDeltaMu2[nadj]);
378: if (ts->vecs_fpu) {
379: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step,ts->vecs_fpu[nadj]);
380: }
381: if (ts->vecs_fpp) {
382: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step,ts->vecs_fpp[nadj]);
383: }
384: }
385: }
386: }
388: if (ts->vecs_sensi2) {
389: VecResetArray(ts->vec_sensip_col);
390: MatDenseRestoreColumn(ts->mat_sensip,&xarr);
391: }
392: th->status = TS_STEP_COMPLETE;
393: return 0;
394: }
396: static PetscErrorCode TSAdjointStep_Theta(TS ts)
397: {
398: TS_Theta *th = (TS_Theta*)ts->data;
399: TS quadts = ts->quadraturets;
400: Vec *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
401: Vec *VecsDeltaLam2 = th->VecsDeltaLam2,*VecsDeltaMu2 = th->VecsDeltaMu2,*VecsSensi2Temp = th->VecsSensi2Temp;
402: PetscInt nadj;
403: Mat J,Jpre,quadJ = NULL,quadJp = NULL;
404: KSP ksp;
405: PetscScalar *xarr;
406: PetscReal adjoint_time_step;
407: PetscReal adjoint_ptime; /* end time of the adjoint time step (ts->ptime is the start time, ususally ts->ptime is larger than adjoint_ptime) */
409: if (th->Theta == 1.) {
410: TSAdjointStepBEuler_Private(ts);
411: return 0;
412: }
413: th->status = TS_STEP_INCOMPLETE;
414: SNESGetKSP(ts->snes,&ksp);
415: TSGetIJacobian(ts,&J,&Jpre,NULL,NULL);
416: if (quadts) {
417: TSGetRHSJacobian(quadts,&quadJ,NULL,NULL,NULL);
418: TSGetRHSJacobianP(quadts,&quadJp,NULL,NULL);
419: }
420: /* If endpoint=1, th->ptime and th->X0 will be used; if endpoint=0, th->stage_time and th->X will be used. */
421: th->stage_time = th->endpoint ? ts->ptime : (ts->ptime+(1.-th->Theta)*ts->time_step);
422: adjoint_ptime = ts->ptime + ts->time_step;
423: adjoint_time_step = -ts->time_step; /* always positive since time_step is negative */
425: if (!th->endpoint) {
426: /* recover th->X0 using vec_sol and the stage value th->X */
427: VecAXPBYPCZ(th->X0,1.0/(1.0-th->Theta),th->Theta/(th->Theta-1.0),0,th->X,ts->vec_sol);
428: }
430: /* Build RHS for first-order adjoint */
431: /* Cost function has an integral term */
432: if (quadts) {
433: if (th->endpoint) {
434: TSComputeRHSJacobian(quadts,th->stage_time,ts->vec_sol,quadJ,NULL);
435: } else {
436: TSComputeRHSJacobian(quadts,th->stage_time,th->X,quadJ,NULL);
437: }
438: }
440: for (nadj=0; nadj<ts->numcost; nadj++) {
441: VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);
442: VecScale(VecsSensiTemp[nadj],1./(th->Theta*adjoint_time_step));
443: if (quadJ) {
444: MatDenseGetColumn(quadJ,nadj,&xarr);
445: VecPlaceArray(ts->vec_drdu_col,xarr);
446: VecAXPY(VecsSensiTemp[nadj],1.,ts->vec_drdu_col);
447: VecResetArray(ts->vec_drdu_col);
448: MatDenseRestoreColumn(quadJ,&xarr);
449: }
450: }
452: /* Build LHS for first-order adjoint */
453: th->shift = 1./(th->Theta*adjoint_time_step);
454: if (th->endpoint) {
455: TSComputeSNESJacobian(ts,ts->vec_sol,J,Jpre);
456: } else {
457: TSComputeSNESJacobian(ts,th->X,J,Jpre);
458: }
459: KSPSetOperators(ksp,J,Jpre);
461: /* Solve stage equation LHS*lambda_s = RHS for first-order adjoint */
462: for (nadj=0; nadj<ts->numcost; nadj++) {
463: KSPConvergedReason kspreason;
464: KSPSolveTranspose(ksp,VecsSensiTemp[nadj],VecsDeltaLam[nadj]);
465: KSPGetConvergedReason(ksp,&kspreason);
466: if (kspreason < 0) {
467: ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
468: PetscInfo(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n",ts->steps,nadj);
469: }
470: }
472: /* Second-order adjoint */
473: if (ts->vecs_sensi2) { /* U_{n+1} */
475: /* Get w1 at t_{n+1} from TLM matrix */
476: MatDenseGetColumn(ts->mat_sensip,0,&xarr);
477: VecPlaceArray(ts->vec_sensip_col,xarr);
478: /* lambda_s^T F_UU w_1 */
479: TSComputeIHessianProductFunctionUU(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fuu);
480: VecResetArray(ts->vec_sensip_col);
481: MatDenseRestoreColumn(ts->mat_sensip,&xarr);
482: /* lambda_s^T F_UP w_2 */
483: TSComputeIHessianProductFunctionUP(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_dir,ts->vecs_fup);
484: for (nadj=0; nadj<ts->numcost; nadj++) { /* compute the residual */
485: VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);
486: VecScale(VecsSensi2Temp[nadj],th->shift);
487: VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fuu[nadj]);
488: if (ts->vecs_fup) {
489: VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fup[nadj]);
490: }
491: }
492: /* Solve stage equation LHS X = RHS for second-order adjoint */
493: for (nadj=0; nadj<ts->numcost; nadj++) {
494: KSPConvergedReason kspreason;
495: KSPSolveTranspose(ksp,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj]);
496: KSPGetConvergedReason(ksp,&kspreason);
497: if (kspreason < 0) {
498: ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
499: PetscInfo(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 2nd-order adjoint solve\n",ts->steps,nadj);
500: }
501: }
502: }
504: /* Update sensitivities, and evaluate integrals if there is any */
505: if (th->endpoint) { /* two-stage Theta methods with th->Theta!=1, th->Theta==1 leads to BEuler */
506: th->shift = 1./((th->Theta-1.)*adjoint_time_step);
507: th->stage_time = adjoint_ptime;
508: TSComputeSNESJacobian(ts,th->X0,J,Jpre);
509: KSPSetOperators(ksp,J,Jpre);
510: /* R_U at t_n */
511: if (quadts) {
512: TSComputeRHSJacobian(quadts,adjoint_ptime,th->X0,quadJ,NULL);
513: }
514: for (nadj=0; nadj<ts->numcost; nadj++) {
515: MatMultTranspose(J,VecsDeltaLam[nadj],ts->vecs_sensi[nadj]);
516: if (quadJ) {
517: MatDenseGetColumn(quadJ,nadj,&xarr);
518: VecPlaceArray(ts->vec_drdu_col,xarr);
519: VecAXPY(ts->vecs_sensi[nadj],-1.,ts->vec_drdu_col);
520: VecResetArray(ts->vec_drdu_col);
521: MatDenseRestoreColumn(quadJ,&xarr);
522: }
523: VecScale(ts->vecs_sensi[nadj],1./th->shift);
524: }
526: /* Second-order adjoint */
527: if (ts->vecs_sensi2) { /* U_n */
528: /* Get w1 at t_n from TLM matrix */
529: MatDenseGetColumn(th->MatFwdSensip0,0,&xarr);
530: VecPlaceArray(ts->vec_sensip_col,xarr);
531: /* lambda_s^T F_UU w_1 */
532: TSComputeIHessianProductFunctionUU(ts,adjoint_ptime,th->X0,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fuu);
533: VecResetArray(ts->vec_sensip_col);
534: MatDenseRestoreColumn(th->MatFwdSensip0,&xarr);
535: /* lambda_s^T F_UU w_2 */
536: TSComputeIHessianProductFunctionUP(ts,adjoint_ptime,th->X0,VecsDeltaLam,ts->vec_dir,ts->vecs_fup);
537: for (nadj=0; nadj<ts->numcost; nadj++) {
538: /* M^T Lambda_s + h(1-theta) F_U^T Lambda_s + h(1-theta) lambda_s^T F_UU w_1 + lambda_s^T F_UP w_2 */
539: MatMultTranspose(J,VecsDeltaLam2[nadj],ts->vecs_sensi2[nadj]);
540: VecAXPY(ts->vecs_sensi2[nadj],1.,ts->vecs_fuu[nadj]);
541: if (ts->vecs_fup) {
542: VecAXPY(ts->vecs_sensi2[nadj],1.,ts->vecs_fup[nadj]);
543: }
544: VecScale(ts->vecs_sensi2[nadj],1./th->shift);
545: }
546: }
548: th->stage_time = ts->ptime; /* recover the old value */
550: if (ts->vecs_sensip) { /* sensitivities wrt parameters */
551: /* U_{n+1} */
552: th->shift = 1.0/(adjoint_time_step*th->Theta);
553: VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->X0,ts->vec_sol);
554: TSComputeIJacobianP(ts,th->stage_time,ts->vec_sol,th->Xdot,-1./(th->Theta*adjoint_time_step),ts->Jacp,PETSC_FALSE);
555: if (quadts) {
556: TSComputeRHSJacobianP(quadts,th->stage_time,ts->vec_sol,quadJp);
557: }
558: for (nadj=0; nadj<ts->numcost; nadj++) {
559: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
560: VecAXPY(ts->vecs_sensip[nadj],-adjoint_time_step*th->Theta,VecsDeltaMu[nadj]);
561: if (quadJp) {
562: MatDenseGetColumn(quadJp,nadj,&xarr);
563: VecPlaceArray(ts->vec_drdp_col,xarr);
564: VecAXPY(ts->vecs_sensip[nadj],adjoint_time_step*th->Theta,ts->vec_drdp_col);
565: VecResetArray(ts->vec_drdp_col);
566: MatDenseRestoreColumn(quadJp,&xarr);
567: }
568: }
569: if (ts->vecs_sensi2p) { /* second-order */
570: /* Get w1 at t_{n+1} from TLM matrix */
571: MatDenseGetColumn(ts->mat_sensip,0,&xarr);
572: VecPlaceArray(ts->vec_sensip_col,xarr);
573: /* lambda_s^T F_PU w_1 */
574: TSComputeIHessianProductFunctionPU(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fpu);
575: VecResetArray(ts->vec_sensip_col);
576: MatDenseRestoreColumn(ts->mat_sensip,&xarr);
578: /* lambda_s^T F_PP w_2 */
579: TSComputeIHessianProductFunctionPP(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_dir,ts->vecs_fpp);
580: for (nadj=0; nadj<ts->numcost; nadj++) {
581: /* Mu2 <- Mu2 + h theta F_P^T Lambda_s + h theta (lambda_s^T F_UU w_1 + lambda_s^T F_UP w_2) */
582: MatMultTranspose(ts->Jacp,VecsDeltaLam2[nadj],VecsDeltaMu2[nadj]);
583: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step*th->Theta,VecsDeltaMu2[nadj]);
584: if (ts->vecs_fpu) {
585: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step*th->Theta,ts->vecs_fpu[nadj]);
586: }
587: if (ts->vecs_fpp) {
588: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step*th->Theta,ts->vecs_fpp[nadj]);
589: }
590: }
591: }
593: /* U_s */
594: VecZeroEntries(th->Xdot);
595: TSComputeIJacobianP(ts,adjoint_ptime,th->X0,th->Xdot,1./((th->Theta-1.0)*adjoint_time_step),ts->Jacp,PETSC_FALSE);
596: if (quadts) {
597: TSComputeRHSJacobianP(quadts,adjoint_ptime,th->X0,quadJp);
598: }
599: for (nadj=0; nadj<ts->numcost; nadj++) {
600: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
601: VecAXPY(ts->vecs_sensip[nadj],-adjoint_time_step*(1.0-th->Theta),VecsDeltaMu[nadj]);
602: if (quadJp) {
603: MatDenseGetColumn(quadJp,nadj,&xarr);
604: VecPlaceArray(ts->vec_drdp_col,xarr);
605: VecAXPY(ts->vecs_sensip[nadj],adjoint_time_step*(1.0-th->Theta),ts->vec_drdp_col);
606: VecResetArray(ts->vec_drdp_col);
607: MatDenseRestoreColumn(quadJp,&xarr);
608: }
609: if (ts->vecs_sensi2p) { /* second-order */
610: /* Get w1 at t_n from TLM matrix */
611: MatDenseGetColumn(th->MatFwdSensip0,0,&xarr);
612: VecPlaceArray(ts->vec_sensip_col,xarr);
613: /* lambda_s^T F_PU w_1 */
614: TSComputeIHessianProductFunctionPU(ts,adjoint_ptime,th->X0,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fpu);
615: VecResetArray(ts->vec_sensip_col);
616: MatDenseRestoreColumn(th->MatFwdSensip0,&xarr);
617: /* lambda_s^T F_PP w_2 */
618: TSComputeIHessianProductFunctionPP(ts,adjoint_ptime,th->X0,VecsDeltaLam,ts->vec_dir,ts->vecs_fpp);
619: for (nadj=0; nadj<ts->numcost; nadj++) {
620: /* Mu2 <- Mu2 + h(1-theta) F_P^T Lambda_s + h(1-theta) (lambda_s^T F_UU w_1 + lambda_s^T F_UP w_2) */
621: MatMultTranspose(ts->Jacp,VecsDeltaLam2[nadj],VecsDeltaMu2[nadj]);
622: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step*(1.0-th->Theta),VecsDeltaMu2[nadj]);
623: if (ts->vecs_fpu) {
624: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step*(1.0-th->Theta),ts->vecs_fpu[nadj]);
625: }
626: if (ts->vecs_fpp) {
627: VecAXPY(ts->vecs_sensi2p[nadj],-adjoint_time_step*(1.0-th->Theta),ts->vecs_fpp[nadj]);
628: }
629: }
630: }
631: }
632: }
633: } else { /* one-stage case */
634: th->shift = 0.0;
635: TSComputeSNESJacobian(ts,th->X,J,Jpre); /* get -f_y */
636: KSPSetOperators(ksp,J,Jpre);
637: if (quadts) {
638: TSComputeRHSJacobian(quadts,th->stage_time,th->X,quadJ,NULL);
639: }
640: for (nadj=0; nadj<ts->numcost; nadj++) {
641: MatMultTranspose(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj]);
642: VecAXPY(ts->vecs_sensi[nadj],-adjoint_time_step,VecsSensiTemp[nadj]);
643: if (quadJ) {
644: MatDenseGetColumn(quadJ,nadj,&xarr);
645: VecPlaceArray(ts->vec_drdu_col,xarr);
646: VecAXPY(ts->vecs_sensi[nadj],adjoint_time_step,ts->vec_drdu_col);
647: VecResetArray(ts->vec_drdu_col);
648: MatDenseRestoreColumn(quadJ,&xarr);
649: }
650: }
651: if (ts->vecs_sensip) {
652: th->shift = 1.0/(adjoint_time_step*th->Theta);
653: VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->X0,th->X);
654: TSComputeIJacobianP(ts,th->stage_time,th->X,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
655: if (quadts) {
656: TSComputeRHSJacobianP(quadts,th->stage_time,th->X,quadJp);
657: }
658: for (nadj=0; nadj<ts->numcost; nadj++) {
659: MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
660: VecAXPY(ts->vecs_sensip[nadj],-adjoint_time_step,VecsDeltaMu[nadj]);
661: if (quadJp) {
662: MatDenseGetColumn(quadJp,nadj,&xarr);
663: VecPlaceArray(ts->vec_drdp_col,xarr);
664: VecAXPY(ts->vecs_sensip[nadj],adjoint_time_step,ts->vec_drdp_col);
665: VecResetArray(ts->vec_drdp_col);
666: MatDenseRestoreColumn(quadJp,&xarr);
667: }
668: }
669: }
670: }
672: th->status = TS_STEP_COMPLETE;
673: return 0;
674: }
676: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
677: {
678: TS_Theta *th = (TS_Theta*)ts->data;
679: PetscReal dt = t - ts->ptime;
681: VecCopy(ts->vec_sol,th->X);
682: if (th->endpoint) dt *= th->Theta;
683: VecWAXPY(X,dt,th->Xdot,th->X);
684: return 0;
685: }
687: static PetscErrorCode TSEvaluateWLTE_Theta(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte)
688: {
689: TS_Theta *th = (TS_Theta*)ts->data;
690: Vec X = ts->vec_sol; /* X = solution */
691: Vec Y = th->vec_lte_work; /* Y = X + LTE */
692: PetscReal wltea,wlter;
694: if (!th->vec_sol_prev) {*wlte = -1; return 0;}
695: /* Cannot compute LTE in first step or in restart after event */
696: if (ts->steprestart) {*wlte = -1; return 0;}
697: /* Compute LTE using backward differences with non-constant time step */
698: {
699: PetscReal h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
700: PetscReal a = 1 + h_prev/h;
701: PetscScalar scal[3]; Vec vecs[3];
702: scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
703: vecs[0] = X; vecs[1] = th->X0; vecs[2] = th->vec_sol_prev;
704: VecCopy(X,Y);
705: VecMAXPY(Y,3,scal,vecs);
706: TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
707: }
708: if (order) *order = 2;
709: return 0;
710: }
712: static PetscErrorCode TSRollBack_Theta(TS ts)
713: {
714: TS_Theta *th = (TS_Theta*)ts->data;
715: TS quadts = ts->quadraturets;
717: VecCopy(th->X0,ts->vec_sol);
718: if (quadts && ts->costintegralfwd) {
719: VecCopy(th->VecCostIntegral0,quadts->vec_sol);
720: }
721: th->status = TS_STEP_INCOMPLETE;
722: if (ts->mat_sensip) {
723: MatCopy(th->MatFwdSensip0,ts->mat_sensip,SAME_NONZERO_PATTERN);
724: }
725: if (quadts && quadts->mat_sensip) {
726: MatCopy(th->MatIntegralSensip0,quadts->mat_sensip,SAME_NONZERO_PATTERN);
727: }
728: return 0;
729: }
731: static PetscErrorCode TSForwardStep_Theta(TS ts)
732: {
733: TS_Theta *th = (TS_Theta*)ts->data;
734: TS quadts = ts->quadraturets;
735: Mat MatDeltaFwdSensip = th->MatDeltaFwdSensip;
736: Vec VecDeltaFwdSensipCol = th->VecDeltaFwdSensipCol;
737: PetscInt ntlm;
738: KSP ksp;
739: Mat J,Jpre,quadJ = NULL,quadJp = NULL;
740: PetscScalar *barr,*xarr;
741: PetscReal previous_shift;
743: previous_shift = th->shift;
744: MatCopy(ts->mat_sensip,th->MatFwdSensip0,SAME_NONZERO_PATTERN);
746: if (quadts && quadts->mat_sensip) {
747: MatCopy(quadts->mat_sensip,th->MatIntegralSensip0,SAME_NONZERO_PATTERN);
748: }
749: SNESGetKSP(ts->snes,&ksp);
750: TSGetIJacobian(ts,&J,&Jpre,NULL,NULL);
751: if (quadts) {
752: TSGetRHSJacobian(quadts,&quadJ,NULL,NULL,NULL);
753: TSGetRHSJacobianP(quadts,&quadJp,NULL,NULL);
754: }
756: /* Build RHS */
757: if (th->endpoint) { /* 2-stage method*/
758: th->shift = 1./((th->Theta-1.)*th->time_step0);
759: TSComputeIJacobian(ts,th->ptime0,th->X0,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
760: MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
761: MatScale(MatDeltaFwdSensip,(th->Theta-1.)/th->Theta);
763: /* Add the f_p forcing terms */
764: if (ts->Jacp) {
765: VecZeroEntries(th->Xdot);
766: TSComputeIJacobianP(ts,th->ptime0,th->X0,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
767: MatAXPY(MatDeltaFwdSensip,(th->Theta-1.)/th->Theta,ts->Jacp,SUBSET_NONZERO_PATTERN);
768: th->shift = previous_shift;
769: VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->X0,ts->vec_sol);
770: TSComputeIJacobianP(ts,th->stage_time,ts->vec_sol,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
771: MatAXPY(MatDeltaFwdSensip,-1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
772: }
773: } else { /* 1-stage method */
774: th->shift = 0.0;
775: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
776: MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
777: MatScale(MatDeltaFwdSensip,-1.);
779: /* Add the f_p forcing terms */
780: if (ts->Jacp) {
781: th->shift = previous_shift;
782: VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->X0,th->X);
783: TSComputeIJacobianP(ts,th->stage_time,th->X,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
784: MatAXPY(MatDeltaFwdSensip,-1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
785: }
786: }
788: /* Build LHS */
789: th->shift = previous_shift; /* recover the previous shift used in TSStep_Theta() */
790: if (th->endpoint) {
791: TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
792: } else {
793: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
794: }
795: KSPSetOperators(ksp,J,Jpre);
797: /*
798: Evaluate the first stage of integral gradients with the 2-stage method:
799: drdu|t_n*S(t_n) + drdp|t_n
800: This is done before the linear solve because the sensitivity variable S(t_n) will be propagated to S(t_{n+1})
801: */
802: if (th->endpoint) { /* 2-stage method only */
803: if (quadts && quadts->mat_sensip) {
804: TSComputeRHSJacobian(quadts,th->ptime0,th->X0,quadJ,NULL);
805: TSComputeRHSJacobianP(quadts,th->ptime0,th->X0,quadJp);
806: MatTransposeMatMult(ts->mat_sensip,quadJ,MAT_REUSE_MATRIX,PETSC_DEFAULT,&th->MatIntegralSensipTemp);
807: MatAXPY(th->MatIntegralSensipTemp,1,quadJp,SAME_NONZERO_PATTERN);
808: MatAXPY(quadts->mat_sensip,th->time_step0*(1.-th->Theta),th->MatIntegralSensipTemp,SAME_NONZERO_PATTERN);
809: }
810: }
812: /* Solve the tangent linear equation for forward sensitivities to parameters */
813: for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
814: KSPConvergedReason kspreason;
815: MatDenseGetColumn(MatDeltaFwdSensip,ntlm,&barr);
816: VecPlaceArray(VecDeltaFwdSensipCol,barr);
817: if (th->endpoint) {
818: MatDenseGetColumn(ts->mat_sensip,ntlm,&xarr);
819: VecPlaceArray(ts->vec_sensip_col,xarr);
820: KSPSolve(ksp,VecDeltaFwdSensipCol,ts->vec_sensip_col);
821: VecResetArray(ts->vec_sensip_col);
822: MatDenseRestoreColumn(ts->mat_sensip,&xarr);
823: } else {
824: KSPSolve(ksp,VecDeltaFwdSensipCol,VecDeltaFwdSensipCol);
825: }
826: KSPGetConvergedReason(ksp,&kspreason);
827: if (kspreason < 0) {
828: ts->reason = TSFORWARD_DIVERGED_LINEAR_SOLVE;
829: PetscInfo(ts,"Step=%D, %Dth tangent linear solve, linear solve fails, stopping tangent linear solve\n",ts->steps,ntlm);
830: }
831: VecResetArray(VecDeltaFwdSensipCol);
832: MatDenseRestoreColumn(MatDeltaFwdSensip,&barr);
833: }
835: /*
836: Evaluate the second stage of integral gradients with the 2-stage method:
837: drdu|t_{n+1}*S(t_{n+1}) + drdp|t_{n+1}
838: */
839: if (quadts && quadts->mat_sensip) {
840: if (!th->endpoint) {
841: MatAXPY(ts->mat_sensip,1,MatDeltaFwdSensip,SAME_NONZERO_PATTERN); /* stage sensitivity */
842: TSComputeRHSJacobian(quadts,th->stage_time,th->X,quadJ,NULL);
843: TSComputeRHSJacobianP(quadts,th->stage_time,th->X,quadJp);
844: MatTransposeMatMult(ts->mat_sensip,quadJ,MAT_REUSE_MATRIX,PETSC_DEFAULT,&th->MatIntegralSensipTemp);
845: MatAXPY(th->MatIntegralSensipTemp,1,quadJp,SAME_NONZERO_PATTERN);
846: MatAXPY(quadts->mat_sensip,th->time_step0,th->MatIntegralSensipTemp,SAME_NONZERO_PATTERN);
847: MatAXPY(ts->mat_sensip,(1.-th->Theta)/th->Theta,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
848: } else {
849: TSComputeRHSJacobian(quadts,th->stage_time,ts->vec_sol,quadJ,NULL);
850: TSComputeRHSJacobianP(quadts,th->stage_time,ts->vec_sol,quadJp);
851: MatTransposeMatMult(ts->mat_sensip,quadJ,MAT_REUSE_MATRIX,PETSC_DEFAULT,&th->MatIntegralSensipTemp);
852: MatAXPY(th->MatIntegralSensipTemp,1,quadJp,SAME_NONZERO_PATTERN);
853: MatAXPY(quadts->mat_sensip,th->time_step0*th->Theta,th->MatIntegralSensipTemp,SAME_NONZERO_PATTERN);
854: }
855: } else {
856: if (!th->endpoint) {
857: MatAXPY(ts->mat_sensip,1./th->Theta,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
858: }
859: }
860: return 0;
861: }
863: static PetscErrorCode TSForwardGetStages_Theta(TS ts,PetscInt *ns,Mat *stagesensip[])
864: {
865: TS_Theta *th = (TS_Theta*)ts->data;
867: if (ns) {
868: if (!th->endpoint && th->Theta != 1.0) *ns = 1; /* midpoint form */
869: else *ns = 2; /* endpoint form */
870: }
871: if (stagesensip) {
872: if (!th->endpoint && th->Theta != 1.0) {
873: th->MatFwdStages[0] = th->MatDeltaFwdSensip;
874: } else {
875: th->MatFwdStages[0] = th->MatFwdSensip0;
876: th->MatFwdStages[1] = ts->mat_sensip; /* stiffly accurate */
877: }
878: *stagesensip = th->MatFwdStages;
879: }
880: return 0;
881: }
883: /*------------------------------------------------------------*/
884: static PetscErrorCode TSReset_Theta(TS ts)
885: {
886: TS_Theta *th = (TS_Theta*)ts->data;
888: VecDestroy(&th->X);
889: VecDestroy(&th->Xdot);
890: VecDestroy(&th->X0);
891: VecDestroy(&th->affine);
893: VecDestroy(&th->vec_sol_prev);
894: VecDestroy(&th->vec_lte_work);
896: VecDestroy(&th->VecCostIntegral0);
897: return 0;
898: }
900: static PetscErrorCode TSAdjointReset_Theta(TS ts)
901: {
902: TS_Theta *th = (TS_Theta*)ts->data;
904: VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
905: VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
906: VecDestroyVecs(ts->numcost,&th->VecsDeltaLam2);
907: VecDestroyVecs(ts->numcost,&th->VecsDeltaMu2);
908: VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);
909: VecDestroyVecs(ts->numcost,&th->VecsSensi2Temp);
910: return 0;
911: }
913: static PetscErrorCode TSDestroy_Theta(TS ts)
914: {
915: TSReset_Theta(ts);
916: if (ts->dm) {
917: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
918: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
919: }
920: PetscFree(ts->data);
921: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
922: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
923: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
924: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
925: return 0;
926: }
928: /*
929: This defines the nonlinear equation that is to be solved with SNES
930: G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
932: Note that U here is the stage argument. This means that U = U_{n+1} only if endpoint = true,
933: otherwise U = theta U_{n+1} + (1 - theta) U0, which for the case of implicit midpoint is
934: U = (U_{n+1} + U0)/2
935: */
936: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
937: {
938: TS_Theta *th = (TS_Theta*)ts->data;
939: Vec X0,Xdot;
940: DM dm,dmsave;
941: PetscReal shift = th->shift;
943: SNESGetDM(snes,&dm);
944: /* When using the endpoint variant, this is actually 1/Theta * Xdot */
945: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
946: if (x != X0) {
947: VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);
948: } else {
949: VecZeroEntries(Xdot);
950: }
951: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
952: dmsave = ts->dm;
953: ts->dm = dm;
954: TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
955: ts->dm = dmsave;
956: TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
957: return 0;
958: }
960: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
961: {
962: TS_Theta *th = (TS_Theta*)ts->data;
963: Vec Xdot;
964: DM dm,dmsave;
965: PetscReal shift = th->shift;
967: SNESGetDM(snes,&dm);
968: /* Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
969: TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);
971: dmsave = ts->dm;
972: ts->dm = dm;
973: TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
974: ts->dm = dmsave;
975: TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
976: return 0;
977: }
979: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
980: {
981: TS_Theta *th = (TS_Theta*)ts->data;
982: TS quadts = ts->quadraturets;
984: /* combine sensitivities to parameters and sensitivities to initial values into one array */
985: th->num_tlm = ts->num_parameters;
986: MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatDeltaFwdSensip);
987: if (quadts && quadts->mat_sensip) {
988: MatDuplicate(quadts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatIntegralSensipTemp);
989: MatDuplicate(quadts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatIntegralSensip0);
990: }
991: /* backup sensitivity results for roll-backs */
992: MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatFwdSensip0);
994: VecDuplicate(ts->vec_sol,&th->VecDeltaFwdSensipCol);
995: return 0;
996: }
998: static PetscErrorCode TSForwardReset_Theta(TS ts)
999: {
1000: TS_Theta *th = (TS_Theta*)ts->data;
1001: TS quadts = ts->quadraturets;
1003: if (quadts && quadts->mat_sensip) {
1004: MatDestroy(&th->MatIntegralSensipTemp);
1005: MatDestroy(&th->MatIntegralSensip0);
1006: }
1007: VecDestroy(&th->VecDeltaFwdSensipCol);
1008: MatDestroy(&th->MatDeltaFwdSensip);
1009: MatDestroy(&th->MatFwdSensip0);
1010: return 0;
1011: }
1013: static PetscErrorCode TSSetUp_Theta(TS ts)
1014: {
1015: TS_Theta *th = (TS_Theta*)ts->data;
1016: TS quadts = ts->quadraturets;
1017: PetscBool match;
1019: if (!th->VecCostIntegral0 && quadts && ts->costintegralfwd) { /* back up cost integral */
1020: VecDuplicate(quadts->vec_sol,&th->VecCostIntegral0);
1021: }
1022: if (!th->X) {
1023: VecDuplicate(ts->vec_sol,&th->X);
1024: }
1025: if (!th->Xdot) {
1026: VecDuplicate(ts->vec_sol,&th->Xdot);
1027: }
1028: if (!th->X0) {
1029: VecDuplicate(ts->vec_sol,&th->X0);
1030: }
1031: if (th->endpoint) {
1032: VecDuplicate(ts->vec_sol,&th->affine);
1033: }
1035: th->order = (th->Theta == 0.5) ? 2 : 1;
1036: th->shift = 1/(th->Theta*ts->time_step);
1038: TSGetDM(ts,&ts->dm);
1039: DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
1040: DMSubDomainHookAdd(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
1042: TSGetAdapt(ts,&ts->adapt);
1043: TSAdaptCandidatesClear(ts->adapt);
1044: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
1045: if (!match) {
1046: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
1047: VecDuplicate(ts->vec_sol,&th->vec_lte_work);
1048: }
1049: TSGetSNES(ts,&ts->snes);
1051: ts->stifflyaccurate = (!th->endpoint && th->Theta != 1.0) ? PETSC_FALSE : PETSC_TRUE;
1052: return 0;
1053: }
1055: /*------------------------------------------------------------*/
1057: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
1058: {
1059: TS_Theta *th = (TS_Theta*)ts->data;
1061: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
1062: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
1063: if (ts->vecs_sensip) {
1064: VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
1065: }
1066: if (ts->vecs_sensi2) {
1067: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam2);
1068: VecDuplicateVecs(ts->vecs_sensi2[0],ts->numcost,&th->VecsSensi2Temp);
1069: /* hack ts to make implicit TS solver work when users provide only explicit versions of callbacks (RHSFunction,RHSJacobian,RHSHessian etc.) */
1070: if (!ts->ihessianproduct_fuu) ts->vecs_fuu = ts->vecs_guu;
1071: if (!ts->ihessianproduct_fup) ts->vecs_fup = ts->vecs_gup;
1072: }
1073: if (ts->vecs_sensi2p) {
1074: VecDuplicateVecs(ts->vecs_sensi2p[0],ts->numcost,&th->VecsDeltaMu2);
1075: /* hack ts to make implicit TS solver work when users provide only explicit versions of callbacks (RHSFunction,RHSJacobian,RHSHessian etc.) */
1076: if (!ts->ihessianproduct_fpu) ts->vecs_fpu = ts->vecs_gpu;
1077: if (!ts->ihessianproduct_fpp) ts->vecs_fpp = ts->vecs_gpp;
1078: }
1079: return 0;
1080: }
1082: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
1083: {
1084: TS_Theta *th = (TS_Theta*)ts->data;
1086: PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
1087: {
1088: PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
1089: PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
1090: PetscOptionsBool("-ts_theta_initial_guess_extrapolate","Extrapolate stage initial guess from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
1091: }
1092: PetscOptionsTail();
1093: return 0;
1094: }
1096: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
1097: {
1098: TS_Theta *th = (TS_Theta*)ts->data;
1099: PetscBool iascii;
1101: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1102: if (iascii) {
1103: PetscViewerASCIIPrintf(viewer," Theta=%g\n",(double)th->Theta);
1104: PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
1105: }
1106: return 0;
1107: }
1109: static PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
1110: {
1111: TS_Theta *th = (TS_Theta*)ts->data;
1113: *theta = th->Theta;
1114: return 0;
1115: }
1117: static PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
1118: {
1119: TS_Theta *th = (TS_Theta*)ts->data;
1122: th->Theta = theta;
1123: th->order = (th->Theta == 0.5) ? 2 : 1;
1124: return 0;
1125: }
1127: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
1128: {
1129: TS_Theta *th = (TS_Theta*)ts->data;
1131: *endpoint = th->endpoint;
1132: return 0;
1133: }
1135: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
1136: {
1137: TS_Theta *th = (TS_Theta*)ts->data;
1139: th->endpoint = flg;
1140: return 0;
1141: }
1143: #if defined(PETSC_HAVE_COMPLEX)
1144: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
1145: {
1146: PetscComplex z = xr + xi*PETSC_i,f;
1147: TS_Theta *th = (TS_Theta*)ts->data;
1148: const PetscReal one = 1.0;
1150: f = (one + (one - th->Theta)*z)/(one - th->Theta*z);
1151: *yr = PetscRealPartComplex(f);
1152: *yi = PetscImaginaryPartComplex(f);
1153: return 0;
1154: }
1155: #endif
1157: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec *Y[])
1158: {
1159: TS_Theta *th = (TS_Theta*)ts->data;
1161: if (ns) {
1162: if (!th->endpoint && th->Theta != 1.0) *ns = 1; /* midpoint form */
1163: else *ns = 2; /* endpoint form */
1164: }
1165: if (Y) {
1166: if (!th->endpoint && th->Theta != 1.0) {
1167: th->Stages[0] = th->X;
1168: } else {
1169: th->Stages[0] = th->X0;
1170: th->Stages[1] = ts->vec_sol; /* stiffly accurate */
1171: }
1172: *Y = th->Stages;
1173: }
1174: return 0;
1175: }
1177: /* ------------------------------------------------------------ */
1178: /*MC
1179: TSTHETA - DAE solver using the implicit Theta method
1181: Level: beginner
1183: Options Database:
1184: + -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
1185: . -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of the Theta method
1186: - -ts_theta_initial_guess_extrapolate <flg> - Extrapolate stage initial guess from previous solution (sometimes unstable)
1188: Notes:
1189: $ -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
1190: $ -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
1191: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)
1193: The endpoint variant of the Theta method and backward Euler can be applied to DAE. The midpoint variant is not suitable for DAEs because it is not stiffly accurate.
1195: The midpoint variant is cast as a 1-stage implicit Runge-Kutta method.
1197: .vb
1198: Theta | Theta
1199: -------------
1200: | 1
1201: .ve
1203: For the default Theta=0.5, this is also known as the implicit midpoint rule.
1205: When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:
1207: .vb
1208: 0 | 0 0
1209: 1 | 1-Theta Theta
1210: -------------------
1211: | 1-Theta Theta
1212: .ve
1214: For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).
1216: To apply a diagonally implicit RK method to DAE, the stage formula
1218: $ Y_i = X + h sum_j a_ij Y'_j
1220: is interpreted as a formula for Y'_i in terms of Y_i and known values (Y'_j, j<i)
1222: .seealso: TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()
1224: M*/
1225: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
1226: {
1227: TS_Theta *th;
1229: ts->ops->reset = TSReset_Theta;
1230: ts->ops->adjointreset = TSAdjointReset_Theta;
1231: ts->ops->destroy = TSDestroy_Theta;
1232: ts->ops->view = TSView_Theta;
1233: ts->ops->setup = TSSetUp_Theta;
1234: ts->ops->adjointsetup = TSAdjointSetUp_Theta;
1235: ts->ops->adjointreset = TSAdjointReset_Theta;
1236: ts->ops->step = TSStep_Theta;
1237: ts->ops->interpolate = TSInterpolate_Theta;
1238: ts->ops->evaluatewlte = TSEvaluateWLTE_Theta;
1239: ts->ops->rollback = TSRollBack_Theta;
1240: ts->ops->setfromoptions = TSSetFromOptions_Theta;
1241: ts->ops->snesfunction = SNESTSFormFunction_Theta;
1242: ts->ops->snesjacobian = SNESTSFormJacobian_Theta;
1243: #if defined(PETSC_HAVE_COMPLEX)
1244: ts->ops->linearstability = TSComputeLinearStability_Theta;
1245: #endif
1246: ts->ops->getstages = TSGetStages_Theta;
1247: ts->ops->adjointstep = TSAdjointStep_Theta;
1248: ts->ops->adjointintegral = TSAdjointCostIntegral_Theta;
1249: ts->ops->forwardintegral = TSForwardCostIntegral_Theta;
1250: ts->default_adapt_type = TSADAPTNONE;
1252: ts->ops->forwardsetup = TSForwardSetUp_Theta;
1253: ts->ops->forwardreset = TSForwardReset_Theta;
1254: ts->ops->forwardstep = TSForwardStep_Theta;
1255: ts->ops->forwardgetstages = TSForwardGetStages_Theta;
1257: ts->usessnes = PETSC_TRUE;
1259: PetscNewLog(ts,&th);
1260: ts->data = (void*)th;
1262: th->VecsDeltaLam = NULL;
1263: th->VecsDeltaMu = NULL;
1264: th->VecsSensiTemp = NULL;
1265: th->VecsSensi2Temp = NULL;
1267: th->extrapolate = PETSC_FALSE;
1268: th->Theta = 0.5;
1269: th->order = 2;
1270: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
1271: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
1272: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
1273: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
1274: return 0;
1275: }
1277: /*@
1278: TSThetaGetTheta - Get the abscissa of the stage in (0,1].
1280: Not Collective
1282: Input Parameter:
1283: . ts - timestepping context
1285: Output Parameter:
1286: . theta - stage abscissa
1288: Note:
1289: Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.
1291: Level: Advanced
1293: .seealso: TSThetaSetTheta()
1294: @*/
1295: PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta)
1296: {
1299: PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
1300: return 0;
1301: }
1303: /*@
1304: TSThetaSetTheta - Set the abscissa of the stage in (0,1].
1306: Not Collective
1308: Input Parameters:
1309: + ts - timestepping context
1310: - theta - stage abscissa
1312: Options Database:
1313: . -ts_theta_theta <theta> - set theta
1315: Level: Intermediate
1317: .seealso: TSThetaGetTheta()
1318: @*/
1319: PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta)
1320: {
1322: PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
1323: return 0;
1324: }
1326: /*@
1327: TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
1329: Not Collective
1331: Input Parameter:
1332: . ts - timestepping context
1334: Output Parameter:
1335: . endpoint - PETSC_TRUE when using the endpoint variant
1337: Level: Advanced
1339: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
1340: @*/
1341: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
1342: {
1345: PetscUseMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
1346: return 0;
1347: }
1349: /*@
1350: TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
1352: Not Collective
1354: Input Parameters:
1355: + ts - timestepping context
1356: - flg - PETSC_TRUE to use the endpoint variant
1358: Options Database:
1359: . -ts_theta_endpoint <flg> - use the endpoint variant
1361: Level: Intermediate
1363: .seealso: TSTHETA, TSCN
1364: @*/
1365: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
1366: {
1368: PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
1369: return 0;
1370: }
1372: /*
1373: * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1374: * The creation functions for these specializations are below.
1375: */
1377: static PetscErrorCode TSSetUp_BEuler(TS ts)
1378: {
1379: TS_Theta *th = (TS_Theta*)ts->data;
1383: TSSetUp_Theta(ts);
1384: return 0;
1385: }
1387: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
1388: {
1389: return 0;
1390: }
1392: /*MC
1393: TSBEULER - ODE solver using the implicit backward Euler method
1395: Level: beginner
1397: Notes:
1398: TSBEULER is equivalent to TSTHETA with Theta=1.0
1400: $ -ts_type theta -ts_theta_theta 1.0
1402: .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA
1404: M*/
1405: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1406: {
1407: TSCreate_Theta(ts);
1408: TSThetaSetTheta(ts,1.0);
1409: TSThetaSetEndpoint(ts,PETSC_FALSE);
1410: ts->ops->setup = TSSetUp_BEuler;
1411: ts->ops->view = TSView_BEuler;
1412: return 0;
1413: }
1415: static PetscErrorCode TSSetUp_CN(TS ts)
1416: {
1417: TS_Theta *th = (TS_Theta*)ts->data;
1421: TSSetUp_Theta(ts);
1422: return 0;
1423: }
1425: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
1426: {
1427: return 0;
1428: }
1430: /*MC
1431: TSCN - ODE solver using the implicit Crank-Nicolson method.
1433: Level: beginner
1435: Notes:
1436: TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.
1438: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint
1440: .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA
1442: M*/
1443: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1444: {
1445: TSCreate_Theta(ts);
1446: TSThetaSetTheta(ts,0.5);
1447: TSThetaSetEndpoint(ts,PETSC_TRUE);
1448: ts->ops->setup = TSSetUp_CN;
1449: ts->ops->view = TSView_CN;
1450: return 0;
1451: }