Actual source code: theta.c
petsc-3.9.4 2018-09-11
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 multiplied with Jacobian transpose */
28: Mat MatDeltaFwdSensip; /* Increment of the forward sensitivity at stage */
29: Vec VecDeltaFwdSensipTemp; /* Working vector for holding one column of the sensitivity matrix */
30: Vec VecDeltaFwdSensipTemp2; /* Additional working vector for endpoint case */
31: Mat MatFwdSensip0; /* backup for roll-backs due to events */
32: Vec VecIntegralSensipTemp; /* Working vector for forward integral sensitivity */
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 TSTheta_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: TSTheta_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 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;
442: if (ts->mat_sensip) {
443: MatCopy(th->MatFwdSensip0,ts->mat_sensip,SAME_NONZERO_PATTERN);
444: }
445: if (ts->vecs_integral_sensip) {
446: for (ncost=0;ncost<ts->numcost;ncost++) {
447: VecCopy(th->VecsIntegralSensip0[ncost],ts->vecs_integral_sensip[ncost]);
448: }
449: }
450: return(0);
451: }
453: static PetscErrorCode TSForwardStep_Theta(TS ts)
454: {
455: TS_Theta *th = (TS_Theta*)ts->data;
456: Mat MatDeltaFwdSensip = th->MatDeltaFwdSensip;
457: Vec VecDeltaFwdSensipTemp = th->VecDeltaFwdSensipTemp,VecDeltaFwdSensipTemp2 = th->VecDeltaFwdSensipTemp2;
458: PetscInt ncost,ntlm;
459: KSP ksp;
460: Mat J,Jp;
461: PetscReal shift;
462: PetscScalar *barr,*xarr;
466: MatCopy(ts->mat_sensip,th->MatFwdSensip0,SAME_NONZERO_PATTERN);
468: for (ncost=0; ncost<ts->numcost; ncost++) {
469: if (ts->vecs_integral_sensip) {
470: VecCopy(ts->vecs_integral_sensip[ncost],th->VecsIntegralSensip0[ncost]);
471: }
472: }
474: SNESGetKSP(ts->snes,&ksp);
475: TSGetIJacobian(ts,&J,&Jp,NULL,NULL);
477: /* Build RHS */
478: if (th->endpoint) { /* 2-stage method*/
479: shift = 1./((th->Theta-1.)*th->time_step);
480: TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,shift,J,Jp,PETSC_FALSE);
481: MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
482: MatScale(MatDeltaFwdSensip,(th->Theta-1.)/th->Theta);
484: /* Add the f_p forcing terms */
485: TSAdjointComputeRHSJacobian(ts,th->ptime,th->X0,ts->Jacp);
486: MatAXPY(MatDeltaFwdSensip,(1.-th->Theta)/th->Theta,ts->Jacp,SUBSET_NONZERO_PATTERN);
487: TSAdjointComputeRHSJacobian(ts,th->stage_time,ts->vec_sol,ts->Jacp);
488: MatAXPY(MatDeltaFwdSensip,1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
489: } else { /* 1-stage method */
490: shift = 0.0;
491: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
492: MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
493: MatScale(MatDeltaFwdSensip,-1.);
495: /* Add the f_p forcing terms */
496: TSAdjointComputeRHSJacobian(ts,th->stage_time,th->X,ts->Jacp);
497: MatAXPY(MatDeltaFwdSensip,1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
498: }
500: /* Build LHS */
501: shift = 1/(th->Theta*th->time_step);
502: if (th->endpoint) {
503: TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,shift,J,Jp,PETSC_FALSE);
504: } else {
505: TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,shift,J,Jp,PETSC_FALSE);
506: }
507: KSPSetOperators(ksp,J,Jp);
509: /*
510: Evaluate the first stage of integral gradients with the 2-stage method:
511: drdy|t_n*S(t_n) + drdp|t_n
512: This is done before the linear solve because the sensitivity variable S(t_n) will be propagated to S(t_{n+1})
513: */
514: if (th->endpoint) { /* 2-stage method only */
515: if (ts->vecs_integral_sensip) {
516: TSAdjointComputeDRDYFunction(ts,th->ptime,th->X0,ts->vecs_drdy);
517: TSAdjointComputeDRDPFunction(ts,th->ptime,th->X0,ts->vecs_drdp);
518: for (ncost=0; ncost<ts->numcost; ncost++) {
519: MatMultTranspose(ts->mat_sensip,ts->vecs_drdy[ncost],th->VecIntegralSensipTemp);
520: VecAXPY(th->VecIntegralSensipTemp,1,ts->vecs_drdp[ncost]);
521: VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*(1.-th->Theta),th->VecIntegralSensipTemp);
522: }
523: }
524: }
526: /* Solve the tangent linear equation for forward sensitivities to parameters */
527: for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
528: MatDenseGetColumn(MatDeltaFwdSensip,ntlm,&barr);
529: VecPlaceArray(VecDeltaFwdSensipTemp,barr);
530: if (th->endpoint) {
531: MatDenseGetColumn(ts->mat_sensip,ntlm,&xarr);
532: VecPlaceArray(VecDeltaFwdSensipTemp2,xarr);
533: KSPSolve(ksp,VecDeltaFwdSensipTemp,VecDeltaFwdSensipTemp2);
534: VecResetArray(VecDeltaFwdSensipTemp2);
535: MatDenseRestoreColumn(ts->mat_sensip,&xarr);
536: } else {
537: KSPSolve(ksp,VecDeltaFwdSensipTemp,VecDeltaFwdSensipTemp);
538: }
539: VecResetArray(VecDeltaFwdSensipTemp);
540: MatDenseRestoreColumn(MatDeltaFwdSensip,&barr);
541: }
543: /*
544: Evaluate the second stage of integral gradients with the 2-stage method:
545: drdy|t_{n+1}*S(t_{n+1}) + drdp|t_{n+1}
546: */
547: if (ts->vecs_integral_sensip) {
548: if (!th->endpoint) {
549: MatAXPY(ts->mat_sensip,1,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
550: TSAdjointComputeDRDYFunction(ts,th->stage_time,th->X,ts->vecs_drdy);
551: TSAdjointComputeDRDPFunction(ts,th->stage_time,th->X,ts->vecs_drdp);
552: for (ncost=0; ncost<ts->numcost; ncost++) {
553: MatMultTranspose(ts->mat_sensip,ts->vecs_drdy[ncost],th->VecIntegralSensipTemp);
554: VecAXPY(th->VecIntegralSensipTemp,1,ts->vecs_drdp[ncost]);
555: VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step,th->VecIntegralSensipTemp);
556: }
557: MatAXPY(ts->mat_sensip,(1.-th->Theta)/th->Theta,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
558: } else {
559: TSAdjointComputeDRDYFunction(ts,th->stage_time,ts->vec_sol,ts->vecs_drdy);
560: TSAdjointComputeDRDPFunction(ts,th->stage_time,ts->vec_sol,ts->vecs_drdp);
561: for (ncost=0; ncost<ts->numcost; ncost++) {
562: MatMultTranspose(ts->mat_sensip,ts->vecs_drdy[ncost],th->VecIntegralSensipTemp);
563: VecAXPY(th->VecIntegralSensipTemp,1,ts->vecs_drdp[ncost]);
564: VecAXPY(ts->vecs_integral_sensip[ncost],th->time_step*th->Theta,th->VecIntegralSensipTemp);
565: }
566: }
567: } else {
568: if (!th->endpoint) {
569: MatAXPY(ts->mat_sensip,1./th->Theta,MatDeltaFwdSensip,SAME_NONZERO_PATTERN);
570: }
571: }
572: return(0);
573: }
575: /*------------------------------------------------------------*/
576: static PetscErrorCode TSReset_Theta(TS ts)
577: {
578: TS_Theta *th = (TS_Theta*)ts->data;
582: VecDestroy(&th->X);
583: VecDestroy(&th->Xdot);
584: VecDestroy(&th->X0);
585: VecDestroy(&th->affine);
587: VecDestroy(&th->vec_sol_prev);
588: VecDestroy(&th->vec_lte_work);
590: VecDestroy(&th->VecCostIntegral0);
591: if (ts->forward_solve) {
592: if (ts->vecs_integral_sensip) {
593: VecDestroy(&th->VecIntegralSensipTemp);
594: VecDestroyVecs(ts->numcost,&th->VecsIntegralSensip0);
595: }
596: VecDestroy(&th->VecDeltaFwdSensipTemp);
597: if (th->endpoint) {
598: VecDestroy(&th->VecDeltaFwdSensipTemp2);
599: }
600: MatDestroy(&th->MatDeltaFwdSensip);
601: MatDestroy(&th->MatFwdSensip0);
602: }
603: VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
604: VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
605: VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);
607: return(0);
608: }
610: static PetscErrorCode TSDestroy_Theta(TS ts)
611: {
615: TSReset_Theta(ts);
616: if (ts->dm) {
617: DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
618: DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
619: }
620: PetscFree(ts->data);
621: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
622: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
623: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
624: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
625: return(0);
626: }
628: /*
629: This defines the nonlinear equation that is to be solved with SNES
630: G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
631: */
632: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
633: {
634: TS_Theta *th = (TS_Theta*)ts->data;
636: Vec X0,Xdot;
637: DM dm,dmsave;
638: PetscReal shift = 1/(th->Theta*ts->time_step);
641: SNESGetDM(snes,&dm);
642: /* When using the endpoint variant, this is actually 1/Theta * Xdot */
643: TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
644: VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);
646: /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
647: dmsave = ts->dm;
648: ts->dm = dm;
649: TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
650: ts->dm = dmsave;
651: TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
652: return(0);
653: }
655: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
656: {
657: TS_Theta *th = (TS_Theta*)ts->data;
659: Vec Xdot;
660: DM dm,dmsave;
661: PetscReal shift = 1/(th->Theta*ts->time_step);
664: SNESGetDM(snes,&dm);
665: /* Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
666: TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);
668: dmsave = ts->dm;
669: ts->dm = dm;
670: TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
671: ts->dm = dmsave;
672: TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
673: return(0);
674: }
676: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
677: {
678: TS_Theta *th = (TS_Theta*)ts->data;
682: /* combine sensitivities to parameters and sensitivities to initial values into one array */
683: th->num_tlm = ts->num_parameters;
684: MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatDeltaFwdSensip);
685: if (ts->vecs_integral_sensip) {
686: VecDuplicate(ts->vecs_integral_sensip[0],&th->VecIntegralSensipTemp);
687: }
688: /* backup sensitivity results for roll-backs */
689: MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatFwdSensip0);
691: if (ts->vecs_integral_sensip) {
692: VecDuplicateVecs(ts->vecs_integral_sensip[0],ts->numcost,&th->VecsIntegralSensip0);
693: }
694: VecDuplicate(ts->vec_sol,&th->VecDeltaFwdSensipTemp);
695: if (th->endpoint) {
696: VecDuplicate(ts->vec_sol,&th->VecDeltaFwdSensipTemp2);
697: }
698: return(0);
699: }
701: static PetscErrorCode TSSetUp_Theta(TS ts)
702: {
703: TS_Theta *th = (TS_Theta*)ts->data;
704: PetscBool match;
708: if (!th->VecCostIntegral0 && ts->vec_costintegral && ts->costintegralfwd) { /* back up cost integral */
709: VecDuplicate(ts->vec_costintegral,&th->VecCostIntegral0);
710: }
711: if (!th->X) {
712: VecDuplicate(ts->vec_sol,&th->X);
713: }
714: if (!th->Xdot) {
715: VecDuplicate(ts->vec_sol,&th->Xdot);
716: }
717: if (!th->X0) {
718: VecDuplicate(ts->vec_sol,&th->X0);
719: }
720: if (th->endpoint) {
721: VecDuplicate(ts->vec_sol,&th->affine);
722: }
724: th->order = (th->Theta == 0.5) ? 2 : 1;
726: TSGetDM(ts,&ts->dm);
727: DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
728: DMSubDomainHookAdd(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
730: TSGetAdapt(ts,&ts->adapt);
731: TSAdaptCandidatesClear(ts->adapt);
732: PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
733: if (!match) {
734: VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
735: VecDuplicate(ts->vec_sol,&th->vec_lte_work);
736: }
737: TSGetSNES(ts,&ts->snes);
738: return(0);
739: }
741: /*------------------------------------------------------------*/
743: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
744: {
745: TS_Theta *th = (TS_Theta*)ts->data;
749: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
750: if(ts->vecs_sensip) {
751: VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
752: }
753: VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
754: return(0);
755: }
757: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
758: {
759: TS_Theta *th = (TS_Theta*)ts->data;
763: PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
764: {
765: PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
766: PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
767: PetscOptionsBool("-ts_theta_initial_guess_extrapolate","Extrapolate stage initial guess from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
768: }
769: PetscOptionsTail();
770: return(0);
771: }
773: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
774: {
775: TS_Theta *th = (TS_Theta*)ts->data;
776: PetscBool iascii;
780: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
781: if (iascii) {
782: PetscViewerASCIIPrintf(viewer," Theta=%g\n",(double)th->Theta);
783: PetscViewerASCIIPrintf(viewer," Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
784: }
785: return(0);
786: }
788: static PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
789: {
790: TS_Theta *th = (TS_Theta*)ts->data;
793: *theta = th->Theta;
794: return(0);
795: }
797: static PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
798: {
799: TS_Theta *th = (TS_Theta*)ts->data;
802: if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
803: th->Theta = theta;
804: th->order = (th->Theta == 0.5) ? 2 : 1;
805: return(0);
806: }
808: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
809: {
810: TS_Theta *th = (TS_Theta*)ts->data;
813: *endpoint = th->endpoint;
814: return(0);
815: }
817: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
818: {
819: TS_Theta *th = (TS_Theta*)ts->data;
822: th->endpoint = flg;
823: return(0);
824: }
826: #if defined(PETSC_HAVE_COMPLEX)
827: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
828: {
829: PetscComplex z = xr + xi*PETSC_i,f;
830: TS_Theta *th = (TS_Theta*)ts->data;
831: const PetscReal one = 1.0;
834: f = (one + (one - th->Theta)*z)/(one - th->Theta*z);
835: *yr = PetscRealPartComplex(f);
836: *yi = PetscImaginaryPartComplex(f);
837: return(0);
838: }
839: #endif
841: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
842: {
843: TS_Theta *th = (TS_Theta*)ts->data;
846: if (ns) *ns = 1;
847: if (Y) *Y = th->endpoint ? &(th->X0) : &(th->X);
848: return(0);
849: }
851: /* ------------------------------------------------------------ */
852: /*MC
853: TSTHETA - DAE solver using the implicit Theta method
855: Level: beginner
857: Options Database:
858: + -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
859: . -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of the Theta method
860: - -ts_theta_initial_guess_extrapolate <flg> - Extrapolate stage initial guess from previous solution (sometimes unstable)
862: Notes:
863: $ -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
864: $ -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
865: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)
867: This method can be applied to DAE.
869: This method is cast as a 1-stage implicit Runge-Kutta method.
871: .vb
872: Theta | Theta
873: -------------
874: | 1
875: .ve
877: For the default Theta=0.5, this is also known as the implicit midpoint rule.
879: When the endpoint variant is chosen, the method becomes a 2-stage method with first stage explicit:
881: .vb
882: 0 | 0 0
883: 1 | 1-Theta Theta
884: -------------------
885: | 1-Theta Theta
886: .ve
888: For the default Theta=0.5, this is the trapezoid rule (also known as Crank-Nicolson, see TSCN).
890: To apply a diagonally implicit RK method to DAE, the stage formula
892: $ Y_i = X + h sum_j a_ij Y'_j
894: is interpreted as a formula for Y'_i in terms of Y_i and known values (Y'_j, j<i)
896: .seealso: TSCreate(), TS, TSSetType(), TSCN, TSBEULER, TSThetaSetTheta(), TSThetaSetEndpoint()
898: M*/
899: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
900: {
901: TS_Theta *th;
905: ts->ops->reset = TSReset_Theta;
906: ts->ops->destroy = TSDestroy_Theta;
907: ts->ops->view = TSView_Theta;
908: ts->ops->setup = TSSetUp_Theta;
909: ts->ops->adjointsetup = TSAdjointSetUp_Theta;
910: ts->ops->step = TSStep_Theta;
911: ts->ops->interpolate = TSInterpolate_Theta;
912: ts->ops->evaluatewlte = TSEvaluateWLTE_Theta;
913: ts->ops->rollback = TSRollBack_Theta;
914: ts->ops->setfromoptions = TSSetFromOptions_Theta;
915: ts->ops->snesfunction = SNESTSFormFunction_Theta;
916: ts->ops->snesjacobian = SNESTSFormJacobian_Theta;
917: #if defined(PETSC_HAVE_COMPLEX)
918: ts->ops->linearstability = TSComputeLinearStability_Theta;
919: #endif
920: ts->ops->getstages = TSGetStages_Theta;
921: ts->ops->adjointstep = TSAdjointStep_Theta;
922: ts->ops->adjointintegral = TSAdjointCostIntegral_Theta;
923: ts->ops->forwardintegral = TSForwardCostIntegral_Theta;
924: ts->default_adapt_type = TSADAPTNONE;
925: ts->ops->forwardsetup = TSForwardSetUp_Theta;
926: ts->ops->forwardstep = TSForwardStep_Theta;
928: ts->usessnes = PETSC_TRUE;
930: PetscNewLog(ts,&th);
931: ts->data = (void*)th;
933: th->VecsDeltaLam = NULL;
934: th->VecsDeltaMu = NULL;
935: th->VecsSensiTemp = NULL;
937: th->extrapolate = PETSC_FALSE;
938: th->Theta = 0.5;
939: th->order = 2;
940: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
941: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
942: PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
943: PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
944: return(0);
945: }
947: /*@
948: TSThetaGetTheta - Get the abscissa of the stage in (0,1].
950: Not Collective
952: Input Parameter:
953: . ts - timestepping context
955: Output Parameter:
956: . theta - stage abscissa
958: Note:
959: Use of this function is normally only required to hack TSTHETA to use a modified integration scheme.
961: Level: Advanced
963: .seealso: TSThetaSetTheta()
964: @*/
965: PetscErrorCode TSThetaGetTheta(TS ts,PetscReal *theta)
966: {
972: PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
973: return(0);
974: }
976: /*@
977: TSThetaSetTheta - Set the abscissa of the stage in (0,1].
979: Not Collective
981: Input Parameter:
982: + ts - timestepping context
983: - theta - stage abscissa
985: Options Database:
986: . -ts_theta_theta <theta>
988: Level: Intermediate
990: .seealso: TSThetaGetTheta()
991: @*/
992: PetscErrorCode TSThetaSetTheta(TS ts,PetscReal theta)
993: {
998: PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
999: return(0);
1000: }
1002: /*@
1003: TSThetaGetEndpoint - Gets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
1005: Not Collective
1007: Input Parameter:
1008: . ts - timestepping context
1010: Output Parameter:
1011: . endpoint - PETSC_TRUE when using the endpoint variant
1013: Level: Advanced
1015: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
1016: @*/
1017: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
1018: {
1024: PetscUseMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
1025: return(0);
1026: }
1028: /*@
1029: TSThetaSetEndpoint - Sets whether to use the endpoint variant of the method (e.g. trapezoid/Crank-Nicolson instead of midpoint rule).
1031: Not Collective
1033: Input Parameter:
1034: + ts - timestepping context
1035: - flg - PETSC_TRUE to use the endpoint variant
1037: Options Database:
1038: . -ts_theta_endpoint <flg>
1040: Level: Intermediate
1042: .seealso: TSTHETA, TSCN
1043: @*/
1044: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
1045: {
1050: PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
1051: return(0);
1052: }
1054: /*
1055: * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1056: * The creation functions for these specializations are below.
1057: */
1059: static PetscErrorCode TSSetUp_BEuler(TS ts)
1060: {
1061: TS_Theta *th = (TS_Theta*)ts->data;
1065: 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");
1066: 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");
1067: TSSetUp_Theta(ts);
1068: return(0);
1069: }
1071: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
1072: {
1074: return(0);
1075: }
1077: /*MC
1078: TSBEULER - ODE solver using the implicit backward Euler method
1080: Level: beginner
1082: Notes:
1083: TSBEULER is equivalent to TSTHETA with Theta=1.0
1085: $ -ts_type theta -ts_theta_theta 1.0
1087: .seealso: TSCreate(), TS, TSSetType(), TSEULER, TSCN, TSTHETA
1089: M*/
1090: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1091: {
1095: TSCreate_Theta(ts);
1096: TSThetaSetTheta(ts,1.0);
1097: TSThetaSetEndpoint(ts,PETSC_FALSE);
1098: ts->ops->setup = TSSetUp_BEuler;
1099: ts->ops->view = TSView_BEuler;
1100: return(0);
1101: }
1103: static PetscErrorCode TSSetUp_CN(TS ts)
1104: {
1105: TS_Theta *th = (TS_Theta*)ts->data;
1109: 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");
1110: 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");
1111: TSSetUp_Theta(ts);
1112: return(0);
1113: }
1115: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
1116: {
1118: return(0);
1119: }
1121: /*MC
1122: TSCN - ODE solver using the implicit Crank-Nicolson method.
1124: Level: beginner
1126: Notes:
1127: TSCN is equivalent to TSTHETA with Theta=0.5 and the "endpoint" option set. I.e.
1129: $ -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint
1131: .seealso: TSCreate(), TS, TSSetType(), TSBEULER, TSTHETA
1133: M*/
1134: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1135: {
1139: TSCreate_Theta(ts);
1140: TSThetaSetTheta(ts,0.5);
1141: TSThetaSetEndpoint(ts,PETSC_TRUE);
1142: ts->ops->setup = TSSetUp_CN;
1143: ts->ops->view = TSView_CN;
1144: return(0);
1145: }