Actual source code: theta.c

petsc-master 2019-10-21
Report Typos and Errors
  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:   PetscReal    shift;
 18:   PetscInt     order;
 19:   PetscBool    endpoint;
 20:   PetscBool    extrapolate;
 21:   TSStepStatus status;
 22:   Vec          VecCostIntegral0;         /* Backup for roll-backs due to events */

 24:   /* context for sensitivity analysis */
 25:   PetscInt     num_tlm;                  /* Total number of tangent linear equations */
 26:   Vec          *VecsDeltaLam;            /* Increment of the adjoint sensitivity w.r.t IC at stage */
 27:   Vec          *VecsDeltaMu;             /* Increment of the adjoint sensitivity w.r.t P at stage */
 28:   Vec          *VecsSensiTemp;           /* Vector to be multiplied with Jacobian transpose */
 29:   Mat          MatDeltaFwdSensip;        /* Increment of the forward sensitivity at stage */
 30:   Vec          VecDeltaFwdSensipCol;     /* Working vector for holding one column of the sensitivity matrix */
 31:   Mat          MatFwdSensip0;            /* backup for roll-backs due to events */
 32:   Mat          MatIntegralSensipTemp;    /* Working vector for forward integral sensitivity */
 33:   Mat          MatIntegralSensip0;       /* backup for roll-backs due to events */
 34:   Vec          *VecsDeltaLam2;           /* Increment of the 2nd-order adjoint sensitivity w.r.t IC at stage */
 35:   Vec          *VecsDeltaMu2;            /* Increment of the 2nd-order adjoint sensitivity w.r.t P at stage */
 36:   Vec          *VecsSensi2Temp;          /* Working vectors that holds the residual for the second-order adjoint */
 37:   Vec          *VecsAffine;              /* Working vectors to store residuals */
 38:   /* context for error estimation */
 39:   Vec          vec_sol_prev;
 40:   Vec          vec_lte_work;
 41: } TS_Theta;

 43: static PetscErrorCode TSThetaGetX0AndXdot(TS ts,DM dm,Vec *X0,Vec *Xdot)
 44: {
 45:   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: {

 67:   if (X0) {
 68:     if (dm && dm != ts->dm) {
 69:       DMRestoreNamedGlobalVector(dm,"TSTheta_X0",X0);
 70:     }
 71:   }
 72:   if (Xdot) {
 73:     if (dm && dm != ts->dm) {
 74:       DMRestoreNamedGlobalVector(dm,"TSTheta_Xdot",Xdot);
 75:     }
 76:   }
 77:   return(0);
 78: }

 80: static PetscErrorCode DMCoarsenHook_TSTheta(DM fine,DM coarse,void *ctx)
 81: {
 83:   return(0);
 84: }

 86: static PetscErrorCode DMRestrictHook_TSTheta(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx)
 87: {
 88:   TS             ts = (TS)ctx;
 90:   Vec            X0,Xdot,X0_c,Xdot_c;

 93:   TSThetaGetX0AndXdot(ts,fine,&X0,&Xdot);
 94:   TSThetaGetX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
 95:   MatRestrict(restrct,X0,X0_c);
 96:   MatRestrict(restrct,Xdot,Xdot_c);
 97:   VecPointwiseMult(X0_c,rscale,X0_c);
 98:   VecPointwiseMult(Xdot_c,rscale,Xdot_c);
 99:   TSThetaRestoreX0AndXdot(ts,fine,&X0,&Xdot);
100:   TSThetaRestoreX0AndXdot(ts,coarse,&X0_c,&Xdot_c);
101:   return(0);
102: }

104: static PetscErrorCode DMSubDomainHook_TSTheta(DM dm,DM subdm,void *ctx)
105: {
107:   return(0);
108: }

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

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

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

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

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

131: static PetscErrorCode TSThetaEvaluateCostIntegral(TS ts)
132: {
133:   TS_Theta       *th = (TS_Theta*)ts->data;
134:   TS             quadts = ts->quadraturets;

138:   if (th->endpoint) {
139:     /* Evolve ts->vec_costintegral to compute integrals */
140:     if (th->Theta!=1.0) {
141:       TSComputeRHSFunction(quadts,th->ptime,th->X0,ts->vec_costintegrand);
142:       VecAXPY(quadts->vec_sol,th->time_step*(1.0-th->Theta),ts->vec_costintegrand);
143:     }
144:     TSComputeRHSFunction(quadts,ts->ptime,ts->vec_sol,ts->vec_costintegrand);
145:     VecAXPY(quadts->vec_sol,th->time_step*th->Theta,ts->vec_costintegrand);
146:   } else {
147:     TSComputeRHSFunction(quadts,th->stage_time,th->X,ts->vec_costintegrand);
148:     VecAXPY(quadts->vec_sol,th->time_step,ts->vec_costintegrand);
149:   }
150:   return(0);
151: }

153: static PetscErrorCode TSForwardCostIntegral_Theta(TS ts)
154: {
155:   TS_Theta       *th = (TS_Theta*)ts->data;
156:   TS             quadts = ts->quadraturets;

160:   /* backup cost integral */
161:   VecCopy(quadts->vec_sol,th->VecCostIntegral0);
162:   TSThetaEvaluateCostIntegral(ts);
163:   return(0);
164: }

166: static PetscErrorCode TSAdjointCostIntegral_Theta(TS ts)
167: {

171:   TSThetaEvaluateCostIntegral(ts);
172:   return(0);
173: }

175: static PetscErrorCode TSTheta_SNESSolve(TS ts,Vec b,Vec x)
176: {
177:   PetscInt       nits,lits;

181:   SNESSolve(ts->snes,b,x);
182:   SNESGetIterationNumber(ts->snes,&nits);
183:   SNESGetLinearSolveIterations(ts->snes,&lits);
184:   ts->snes_its += nits; ts->ksp_its += lits;
185:   return(0);
186: }

188: static PetscErrorCode TSStep_Theta(TS ts)
189: {
190:   TS_Theta       *th = (TS_Theta*)ts->data;
191:   PetscInt       rejections = 0;
192:   PetscBool      stageok,accept = PETSC_TRUE;
193:   PetscReal      next_time_step = ts->time_step;

197:   if (!ts->steprollback) {
198:     if (th->vec_sol_prev) { VecCopy(th->X0,th->vec_sol_prev); }
199:     VecCopy(ts->vec_sol,th->X0);
200:   }

202:   th->status = TS_STEP_INCOMPLETE;
203:   while (!ts->reason && th->status != TS_STEP_COMPLETE) {

205:     th->shift = 1/(th->Theta*ts->time_step);
206:     th->stage_time = ts->ptime + (th->endpoint ? (PetscReal)1 : th->Theta)*ts->time_step;

208:     VecCopy(th->X0,th->X);
209:     if (th->extrapolate && !ts->steprestart) {
210:       VecAXPY(th->X,1/th->shift,th->Xdot);
211:     }
212:     if (th->endpoint) { /* This formulation assumes linear time-independent mass matrix */
213:       if (!th->affine) {VecDuplicate(ts->vec_sol,&th->affine);}
214:       VecZeroEntries(th->Xdot);
215:       TSComputeIFunction(ts,ts->ptime,th->X0,th->Xdot,th->affine,PETSC_FALSE);
216:       VecScale(th->affine,(th->Theta-1)/th->Theta);
217:     } else if (th->affine) { /* Just in case th->endpoint is changed between calls to TSStep_Theta() */
218:       VecZeroEntries(th->affine);
219:     }
220:     TSPreStage(ts,th->stage_time);
221:     TSTheta_SNESSolve(ts,th->affine,th->X);
222:     TSPostStage(ts,th->stage_time,0,&th->X);
223:     TSAdaptCheckStage(ts->adapt,ts,th->stage_time,th->X,&stageok);
224:     if (!stageok) goto reject_step;

226:     th->status = TS_STEP_PENDING;
227:     if (th->endpoint) {
228:       VecCopy(th->X,ts->vec_sol);
229:     } else {
230:       VecAXPBYPCZ(th->Xdot,-th->shift,th->shift,0,th->X0,th->X);
231:       VecAXPY(ts->vec_sol,ts->time_step,th->Xdot);
232:     }
233:     TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);
234:     th->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE;
235:     if (!accept) {
236:       VecCopy(th->X0,ts->vec_sol);
237:       ts->time_step = next_time_step;
238:       goto reject_step;
239:     }

241:     if (ts->forward_solve || ts->costintegralfwd) { /* Save the info for the later use in cost integral evaluation */
242:       th->ptime     = ts->ptime;
243:       th->time_step = ts->time_step;
244:     }

246:     ts->ptime += ts->time_step;
247:     ts->time_step = next_time_step;
248:     break;

250:   reject_step:
251:     ts->reject++; accept = PETSC_FALSE;
252:     if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) {
253:       ts->reason = TS_DIVERGED_STEP_REJECTED;
254:       PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);
255:     }
256:   }
257:   return(0);
258: }

260: static PetscErrorCode TSAdjointStepBEuler_Private(TS ts)
261: {
262:   TS_Theta       *th = (TS_Theta*)ts->data;
263:   TS             quadts = ts->quadraturets;
264:   Vec            *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
265:   Vec            *VecsDeltaLam2 = th->VecsDeltaLam2,*VecsDeltaMu2 = th->VecsDeltaMu2,*VecsSensi2Temp = th->VecsSensi2Temp;
266:   PetscInt       nadj;
267:   Mat            J,Jpre,quadJ = NULL,quadJp = NULL;
268:   KSP            ksp;
269:   PetscScalar    *xarr;
270:   TSEquationType eqtype;
271:   PetscBool      isexplicitode = PETSC_FALSE;

275:   TSGetEquationType(ts,&eqtype);
276:   if (eqtype == TS_EQ_ODE_EXPLICIT) {
277:     isexplicitode  = PETSC_TRUE;
278:     VecsDeltaLam  = ts->vecs_sensi;
279:     VecsDeltaLam2 = ts->vecs_sensi2;
280:   }
281:   th->status = TS_STEP_INCOMPLETE;
282:   SNESGetKSP(ts->snes,&ksp);
283:   TSGetIJacobian(ts,&J,&Jpre,NULL,NULL);
284:   if (quadts) {
285:     TSGetRHSJacobian(quadts,&quadJ,NULL,NULL,NULL);
286:     TSGetRHSJacobianP(quadts,&quadJp,NULL,NULL);
287:   }

289:   /* If endpoint=1, th->ptime and th->X0 will be used; if endpoint=0, th->stage_time and th->X will be used. */
290:   th->stage_time = ts->ptime; /* time_step is negative*/
291:   th->ptime      = ts->ptime + ts->time_step;
292:   th->time_step  = -ts->time_step;

294:   /* Build RHS for first-order adjoint lambda_{n+1}/h + r_u^T(n+1) */
295:   if (quadts) {
296:     TSComputeRHSJacobian(quadts,th->stage_time,th->X,quadJ,NULL);
297:   }

299:   for (nadj=0; nadj<ts->numcost; nadj++) {
300:     VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);
301:     VecScale(VecsSensiTemp[nadj],1./th->time_step); /* lambda_{n+1}/h */
302:     if (quadJ) {
303:       MatDenseGetColumn(quadJ,nadj,&xarr);
304:       VecPlaceArray(ts->vec_drdu_col,xarr);
305:       VecAXPY(VecsSensiTemp[nadj],1.,ts->vec_drdu_col);
306:       VecResetArray(ts->vec_drdu_col);
307:       MatDenseRestoreColumn(quadJ,&xarr);
308:     }
309:   }

311:   /* Build LHS for first-order adjoint */
312:   TSComputeSNESJacobian(ts,th->X,J,Jpre);
313:   KSPSetOperators(ksp,J,Jpre);

315:   /* Solve stage equation LHS*lambda_s = RHS for first-order adjoint */
316:   for (nadj=0; nadj<ts->numcost; nadj++) {
317:     KSPConvergedReason kspreason;
318:     KSPSolveTranspose(ksp,VecsSensiTemp[nadj],VecsDeltaLam[nadj]);
319:     KSPGetConvergedReason(ksp,&kspreason);
320:     if (kspreason < 0) {
321:       ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
322:       PetscInfo2(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n",ts->steps,nadj);
323:     }
324:   }

326:   if (ts->vecs_sensi2) { /* U_{n+1} */
327:     /* Get w1 at t_{n+1} from TLM matrix */
328:     MatDenseGetColumn(ts->mat_sensip,0,&xarr);
329:     VecPlaceArray(ts->vec_sensip_col,xarr);
330:     /* lambda_s^T F_UU w_1 */
331:     TSComputeIHessianProductFunctionUU(ts,th->stage_time,th->X,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fuu);
332:     /* lambda_s^T F_UP w_2 */
333:     TSComputeIHessianProductFunctionUP(ts,th->stage_time,th->X,VecsDeltaLam,ts->vec_dir,ts->vecs_fup);
334:     for (nadj=0; nadj<ts->numcost; nadj++) { /* compute the residual */
335:       VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);
336:       VecScale(VecsSensi2Temp[nadj],1./th->time_step);
337:       VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fuu[nadj]);
338:       if (ts->vecs_fup) {
339:         VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fup[nadj]);
340:       }
341:     }
342:     /* Solve stage equation LHS X = RHS for second-order adjoint */
343:     for (nadj=0; nadj<ts->numcost; nadj++) {
344:       KSPConvergedReason kspreason;
345:       KSPSolveTranspose(ksp,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj]);
346:       KSPGetConvergedReason(ksp,&kspreason);
347:       if (kspreason < 0) {
348:         ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
349:         PetscInfo2(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 2nd-order adjoint solve\n",ts->steps,nadj);
350:       }
351:     }
352:   }

354:   /* Update sensitivities, and evaluate integrals if there is any */
355:   if (!isexplicitode) {
356:     th->shift = 0.0;
357:     TSComputeSNESJacobian(ts,th->X,J,Jpre);
358:     KSPSetOperators(ksp,J,Jpre);
359:     MatScale(J,-1.);
360:     for (nadj=0; nadj<ts->numcost; nadj++) {
361:       /* Add f_U \lambda_s to the original RHS */
362:       MatMultTransposeAdd(J,VecsDeltaLam[nadj],VecsSensiTemp[nadj],VecsSensiTemp[nadj]);
363:       VecScale(VecsSensiTemp[nadj],th->time_step);
364:       VecCopy(VecsSensiTemp[nadj],ts->vecs_sensi[nadj]);
365:       if (ts->vecs_sensi2) {
366:         MatMultTransposeAdd(J,VecsDeltaLam2[nadj],VecsSensi2Temp[nadj],VecsSensi2Temp[nadj]);
367:         VecScale(VecsSensi2Temp[nadj],th->time_step);
368:         VecCopy(VecsSensi2Temp[nadj],ts->vecs_sensi2[nadj]);
369:       }
370:     }
371:   }
372:   if (ts->vecs_sensip) {
373:     th->shift = 1./th->time_step;
374:     TSComputeIJacobianP(ts,th->stage_time,th->X,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE); /* get -f_p */
375:     if (quadts) {
376:       TSComputeRHSJacobianP(quadts,th->stage_time,th->X,quadJp);
377:     }
378:     if (ts->vecs_sensi2p) {
379:       /* lambda_s^T F_PU w_1 */
380:       TSComputeIHessianProductFunctionPU(ts,th->stage_time,th->X,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fpu);
381:       /* lambda_s^T F_PP w_2 */
382:       TSComputeIHessianProductFunctionPP(ts,th->stage_time,th->X,VecsDeltaLam,ts->vec_dir,ts->vecs_fpp);
383:     }

385:     for (nadj=0; nadj<ts->numcost; nadj++) {
386:       MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
387:       VecAXPY(ts->vecs_sensip[nadj],-th->time_step,VecsDeltaMu[nadj]);
388:       if (quadJp) {
389:         MatDenseGetColumn(quadJp,nadj,&xarr);
390:         VecPlaceArray(ts->vec_drdp_col,xarr);
391:         VecAXPY(ts->vecs_sensip[nadj],th->time_step,ts->vec_drdp_col);
392:         VecResetArray(ts->vec_drdp_col);
393:         MatDenseRestoreColumn(quadJp,&xarr);
394:       }
395:       if (ts->vecs_sensi2p) {
396:         MatMultTranspose(ts->Jacp,VecsDeltaLam2[nadj],VecsDeltaMu2[nadj]);
397:         VecAXPY(ts->vecs_sensi2p[nadj],-th->time_step,VecsDeltaMu2[nadj]);
398:         if (ts->vecs_fpu) {
399:           VecAXPY(ts->vecs_sensi2p[nadj],-th->time_step,ts->vecs_fpu[nadj]);
400:         }
401:         if (ts->vecs_fpp) {
402:           VecAXPY(ts->vecs_sensi2p[nadj],-th->time_step,ts->vecs_fpp[nadj]);
403:         }
404:       }
405:     }
406:   }

408:   if (ts->vecs_sensi2) {
409:     VecResetArray(ts->vec_sensip_col);
410:     MatDenseRestoreColumn(ts->mat_sensip,&xarr);
411:   }
412:   th->status = TS_STEP_COMPLETE;
413:   return(0);
414: }

416: static PetscErrorCode TSAdjointStep_Theta(TS ts)
417: {
418:   TS_Theta       *th = (TS_Theta*)ts->data;
419:   TS             quadts = ts->quadraturets;
420:   Vec            *VecsDeltaLam = th->VecsDeltaLam,*VecsDeltaMu = th->VecsDeltaMu,*VecsSensiTemp = th->VecsSensiTemp;
421:   Vec            *VecsDeltaLam2 = th->VecsDeltaLam2,*VecsDeltaMu2 = th->VecsDeltaMu2,*VecsSensi2Temp = th->VecsSensi2Temp;
422:   PetscInt       nadj;
423:   Mat            J,Jpre,quadJ = NULL,quadJp = NULL;
424:   KSP            ksp;
425:   PetscScalar    *xarr;

429:   if (th->Theta == 1.) {
430:     TSAdjointStepBEuler_Private(ts);
431:     return(0);
432:   }
433:   th->status = TS_STEP_INCOMPLETE;
434:   SNESGetKSP(ts->snes,&ksp);
435:   TSGetIJacobian(ts,&J,&Jpre,NULL,NULL);
436:   if (quadts) {
437:     TSGetRHSJacobian(quadts,&quadJ,NULL,NULL,NULL);
438:     TSGetRHSJacobianP(quadts,&quadJp,NULL,NULL);
439:   }
440:   /* If endpoint=1, th->ptime and th->X0 will be used; if endpoint=0, th->stage_time and th->X will be used. */
441:   th->stage_time = th->endpoint ? ts->ptime : (ts->ptime+(1.-th->Theta)*ts->time_step); /* time_step is negative*/
442:   th->ptime      = ts->ptime + ts->time_step;
443:   th->time_step  = -ts->time_step;

445:   /* Build RHS for first-order adjoint */
446:   /* Cost function has an integral term */
447:   if (quadts) {
448:     if (th->endpoint) {
449:       TSComputeRHSJacobian(quadts,th->stage_time,ts->vec_sol,quadJ,NULL);
450:     } else {
451:       TSComputeRHSJacobian(quadts,th->stage_time,th->X,quadJ,NULL);
452:     }
453:   }

455:   for (nadj=0; nadj<ts->numcost; nadj++) {
456:     VecCopy(ts->vecs_sensi[nadj],VecsSensiTemp[nadj]);
457:     VecScale(VecsSensiTemp[nadj],1./(th->Theta*th->time_step));
458:     if (quadJ) {
459:       MatDenseGetColumn(quadJ,nadj,&xarr);
460:       VecPlaceArray(ts->vec_drdu_col,xarr);
461:       VecAXPY(VecsSensiTemp[nadj],1.,ts->vec_drdu_col);
462:       VecResetArray(ts->vec_drdu_col);
463:       MatDenseRestoreColumn(quadJ,&xarr);
464:     }
465:   }

467:   /* Build LHS for first-order adjoint */
468:   th->shift = 1./(th->Theta*th->time_step);
469:   if (th->endpoint) {
470:     TSComputeSNESJacobian(ts,ts->vec_sol,J,Jpre);
471:   } else {
472:     TSComputeSNESJacobian(ts,th->X,J,Jpre);
473:   }
474:   KSPSetOperators(ksp,J,Jpre);

476:   /* Solve stage equation LHS*lambda_s = RHS for first-order adjoint */
477:   for (nadj=0; nadj<ts->numcost; nadj++) {
478:     KSPConvergedReason kspreason;
479:     KSPSolveTranspose(ksp,VecsSensiTemp[nadj],VecsDeltaLam[nadj]);
480:     KSPGetConvergedReason(ksp,&kspreason);
481:     if (kspreason < 0) {
482:       ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
483:       PetscInfo2(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 1st-order adjoint solve\n",ts->steps,nadj);
484:     }
485:   }

487:   /* Second-order adjoint */
488:   if (ts->vecs_sensi2) { /* U_{n+1} */
489:     if (!th->endpoint) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Operation not implemented in TS_Theta");
490:     /* Get w1 at t_{n+1} from TLM matrix */
491:     MatDenseGetColumn(ts->mat_sensip,0,&xarr);
492:     VecPlaceArray(ts->vec_sensip_col,xarr);
493:     /* lambda_s^T F_UU w_1 */
494:     TSComputeIHessianProductFunctionUU(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fuu);
495:     VecResetArray(ts->vec_sensip_col);
496:     MatDenseRestoreColumn(ts->mat_sensip,&xarr);
497:     /* lambda_s^T F_UP w_2 */
498:     TSComputeIHessianProductFunctionUP(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_dir,ts->vecs_fup);
499:     for (nadj=0; nadj<ts->numcost; nadj++) { /* compute the residual */
500:       VecCopy(ts->vecs_sensi2[nadj],VecsSensi2Temp[nadj]);
501:       VecScale(VecsSensi2Temp[nadj],th->shift);
502:       VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fuu[nadj]);
503:       if (ts->vecs_fup) {
504:         VecAXPY(VecsSensi2Temp[nadj],-1.,ts->vecs_fup[nadj]);
505:       }
506:     }
507:     /* Solve stage equation LHS X = RHS for second-order adjoint */
508:     for (nadj=0; nadj<ts->numcost; nadj++) {
509:       KSPConvergedReason kspreason;
510:       KSPSolveTranspose(ksp,VecsSensi2Temp[nadj],VecsDeltaLam2[nadj]);
511:       KSPGetConvergedReason(ksp,&kspreason);
512:       if (kspreason < 0) {
513:         ts->reason = TSADJOINT_DIVERGED_LINEAR_SOLVE;
514:         PetscInfo2(ts,"Step=%D, %Dth cost function, transposed linear solve fails, stopping 2nd-order adjoint solve\n",ts->steps,nadj);
515:       }
516:     }
517:   }

519:   /* Update sensitivities, and evaluate integrals if there is any */
520:   if(th->endpoint) { /* two-stage Theta methods with th->Theta!=1, th->Theta==1 leads to BEuler */
521:     th->shift      = 1./((th->Theta-1.)*th->time_step);
522:     th->stage_time = th->ptime;
523:     TSComputeSNESJacobian(ts,th->X0,J,Jpre);
524:     KSPSetOperators(ksp,J,Jpre);
525:     /* R_U at t_n */
526:     if (quadts) {
527:       TSComputeRHSJacobian(quadts,th->ptime,th->X0,quadJ,NULL);
528:     }
529:     for (nadj=0; nadj<ts->numcost; nadj++) {
530:       MatMultTranspose(J,VecsDeltaLam[nadj],ts->vecs_sensi[nadj]);
531:       if (quadJ) {
532:         MatDenseGetColumn(quadJ,nadj,&xarr);
533:         VecPlaceArray(ts->vec_drdu_col,xarr);
534:         VecAXPY(ts->vecs_sensi[nadj],-1.,ts->vec_drdu_col);
535:         VecResetArray(ts->vec_drdu_col);
536:         MatDenseRestoreColumn(quadJ,&xarr);
537:       }
538:       VecScale(ts->vecs_sensi[nadj],1./th->shift);
539:     }

541:     /* Second-order adjoint */
542:     if (ts->vecs_sensi2) { /* U_n */
543:       /* Get w1 at t_n from TLM matrix */
544:       MatDenseGetColumn(th->MatFwdSensip0,0,&xarr);
545:       VecPlaceArray(ts->vec_sensip_col,xarr);
546:       /* lambda_s^T F_UU w_1 */
547:       TSComputeIHessianProductFunctionUU(ts,th->ptime,th->X0,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fuu);
548:       VecResetArray(ts->vec_sensip_col);
549:       MatDenseRestoreColumn(th->MatFwdSensip0,&xarr);
550:       /* lambda_s^T F_UU w_2 */
551:       TSComputeIHessianProductFunctionUP(ts,th->ptime,th->X0,VecsDeltaLam,ts->vec_dir,ts->vecs_fup);
552:       for (nadj=0; nadj<ts->numcost; nadj++) {
553:         /* 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  */
554:         MatMultTranspose(J,VecsDeltaLam2[nadj],ts->vecs_sensi2[nadj]);
555:         VecAXPY(ts->vecs_sensi2[nadj],1.,ts->vecs_fuu[nadj]);
556:         if (ts->vecs_fup) {
557:           VecAXPY(ts->vecs_sensi2[nadj],1.,ts->vecs_fup[nadj]);
558:         }
559:         VecScale(ts->vecs_sensi2[nadj],1./th->shift);
560:       }
561:     }

563:     th->stage_time = ts->ptime; /* recover the old value */

565:     if (ts->vecs_sensip) { /* sensitivities wrt parameters */
566:       /* U_{n+1} */
567:       th->shift = -1./(th->Theta*th->time_step);
568:       TSComputeIJacobianP(ts,th->stage_time,ts->vec_sol,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
569:       if (quadts) {
570:         TSComputeRHSJacobianP(quadts,th->stage_time,ts->vec_sol,quadJp);
571:       }
572:       for (nadj=0; nadj<ts->numcost; nadj++) {
573:         MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
574:         VecAXPY(ts->vecs_sensip[nadj],-th->time_step*th->Theta,VecsDeltaMu[nadj]);
575:       }
576:       if (ts->vecs_sensi2p) { /* second-order */
577:         /* Get w1 at t_{n+1} from TLM matrix */
578:         MatDenseGetColumn(ts->mat_sensip,0,&xarr);
579:         VecPlaceArray(ts->vec_sensip_col,xarr);
580:         /* lambda_s^T F_PU w_1 */
581:         TSComputeIHessianProductFunctionPU(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_sensip_col,ts->vecs_fpu);
582:         VecResetArray(ts->vec_sensip_col);
583:         MatDenseRestoreColumn(ts->mat_sensip,&xarr);

585:         /* lambda_s^T F_PP w_2 */
586:         TSComputeIHessianProductFunctionPP(ts,th->stage_time,ts->vec_sol,VecsDeltaLam,ts->vec_dir,ts->vecs_fpp);
587:         for (nadj=0; nadj<ts->numcost; nadj++) {
588:           /* 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)  */
589:           MatMultTranspose(ts->Jacp,VecsDeltaLam2[nadj],VecsDeltaMu2[nadj]);
590:           VecAXPY(ts->vecs_sensi2p[nadj],-th->time_step*th->Theta,VecsDeltaMu2[nadj]);
591:           if (ts->vecs_fpu) {
592:             VecAXPY(ts->vecs_sensi2p[nadj],-th->time_step*th->Theta,ts->vecs_fpu[nadj]);
593:           }
594:           if (ts->vecs_fpp) {
595:             VecAXPY(ts->vecs_sensi2p[nadj],-th->time_step*th->Theta,ts->vecs_fpp[nadj]);
596:           }
597:         }
598:       }

600:       /* U_s */
601:       th->shift = 1./((th->Theta-1.0)*th->time_step);
602:       TSComputeIJacobianP(ts,th->ptime,th->X0,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
603:       if (quadts) {
604:         TSComputeRHSJacobianP(quadts,th->ptime,th->X0,quadJp);
605:       }
606:       for (nadj=0; nadj<ts->numcost; nadj++) {
607:         MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
608:         VecAXPY(ts->vecs_sensip[nadj],-th->time_step*(1.0-th->Theta),VecsDeltaMu[nadj]);
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,th->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,th->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],-th->time_step*(1.0-th->Theta),VecsDeltaMu2[nadj]);
623:             if (ts->vecs_fpu) {
624:               VecAXPY(ts->vecs_sensi2p[nadj],-th->time_step*(1.0-th->Theta),ts->vecs_fpu[nadj]);
625:             }
626:             if (ts->vecs_fpp) {
627:               VecAXPY(ts->vecs_sensi2p[nadj],-th->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],-th->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],th->time_step,ts->vec_drdu_col);
647:         VecResetArray(ts->vec_drdu_col);
648:         MatDenseRestoreColumn(quadJ,&xarr);
649:       }
650:     }
651:     if (ts->vecs_sensip) {
652:       TSComputeIJacobianP(ts,th->stage_time,th->X,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
653:       if (quadts) {
654:         TSComputeRHSJacobianP(quadts,th->stage_time,th->X,quadJp);
655:       }
656:       for (nadj=0; nadj<ts->numcost; nadj++) {
657:         MatMultTranspose(ts->Jacp,VecsDeltaLam[nadj],VecsDeltaMu[nadj]);
658:         VecAXPY(ts->vecs_sensip[nadj],-th->time_step,VecsDeltaMu[nadj]);
659:         if (quadJp) {
660:           MatDenseGetColumn(quadJp,nadj,&xarr);
661:           VecPlaceArray(ts->vec_drdp_col,xarr);
662:           VecAXPY(ts->vecs_sensip[nadj],th->time_step,ts->vec_drdp_col);
663:           VecResetArray(ts->vec_drdp_col);
664:           MatDenseRestoreColumn(quadJp,&xarr);
665:         }
666:       }
667:     }
668:   }

670:   th->status = TS_STEP_COMPLETE;
671:   return(0);
672: }

674: static PetscErrorCode TSInterpolate_Theta(TS ts,PetscReal t,Vec X)
675: {
676:   TS_Theta       *th = (TS_Theta*)ts->data;
677:   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;

696:   if (!th->vec_sol_prev) {*wlte = -1; return(0);}
697:   /* Cannot compute LTE in first step or in restart after event */
698:   if (ts->steprestart) {*wlte = -1; return(0);}
699:   /* Compute LTE using backward differences with non-constant time step */
700:   {
701:     PetscReal   h = ts->time_step, h_prev = ts->ptime - ts->ptime_prev;
702:     PetscReal   a = 1 + h_prev/h;
703:     PetscScalar scal[3]; Vec vecs[3];
704:     scal[0] = +1/a; scal[1] = -1/(a-1); scal[2] = +1/(a*(a-1));
705:     vecs[0] = X;    vecs[1] = th->X0;   vecs[2] = th->vec_sol_prev;
706:     VecCopy(X,Y);
707:     VecMAXPY(Y,3,scal,vecs);
708:     TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);
709:   }
710:   if (order) *order = 2;
711:   return(0);
712: }

714: static PetscErrorCode TSRollBack_Theta(TS ts)
715: {
716:   TS_Theta       *th = (TS_Theta*)ts->data;
717:   TS             quadts = ts->quadraturets;

721:   VecCopy(th->X0,ts->vec_sol);
722:   if (quadts && ts->costintegralfwd) {
723:     VecCopy(th->VecCostIntegral0,quadts->vec_sol);
724:   }
725:   th->status = TS_STEP_INCOMPLETE;
726:   if (ts->mat_sensip) {
727:     MatCopy(th->MatFwdSensip0,ts->mat_sensip,SAME_NONZERO_PATTERN);
728:   }
729:   if (quadts && quadts->mat_sensip) {
730:     MatCopy(th->MatIntegralSensip0,quadts->mat_sensip,SAME_NONZERO_PATTERN);
731:   }
732:   return(0);
733: }

735: static PetscErrorCode TSForwardStep_Theta(TS ts)
736: {
737:   TS_Theta       *th = (TS_Theta*)ts->data;
738:   TS             quadts = ts->quadraturets;
739:   Mat            MatDeltaFwdSensip = th->MatDeltaFwdSensip;
740:   Vec            VecDeltaFwdSensipCol = th->VecDeltaFwdSensipCol;
741:   PetscInt       ntlm;
742:   KSP            ksp;
743:   Mat            J,Jpre,quadJ = NULL,quadJp = NULL;
744:   PetscScalar    *barr,*xarr;

748:   MatCopy(ts->mat_sensip,th->MatFwdSensip0,SAME_NONZERO_PATTERN);

750:   if (quadts && quadts->mat_sensip) {
751:     MatCopy(quadts->mat_sensip,th->MatIntegralSensip0,SAME_NONZERO_PATTERN);
752:   }
753:   SNESGetKSP(ts->snes,&ksp);
754:   TSGetIJacobian(ts,&J,&Jpre,NULL,NULL);
755:   if (quadts) {
756:     TSGetRHSJacobian(quadts,&quadJ,NULL,NULL,NULL);
757:     TSGetRHSJacobianP(quadts,&quadJp,NULL,NULL);
758:   }

760:   /* Build RHS */
761:   if (th->endpoint) { /* 2-stage method*/
762:     th->shift = 1./((th->Theta-1.)*th->time_step);
763:     TSComputeIJacobian(ts,th->ptime,th->X0,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
764:     MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
765:     MatScale(MatDeltaFwdSensip,(th->Theta-1.)/th->Theta);

767:     /* Add the f_p forcing terms */
768:     if (ts->Jacp) {
769:       TSComputeIJacobianP(ts,th->ptime,th->X0,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
770:       MatAXPY(MatDeltaFwdSensip,(th->Theta-1.)/th->Theta,ts->Jacp,SUBSET_NONZERO_PATTERN);
771:       TSComputeIJacobianP(ts,th->stage_time,ts->vec_sol,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
772:       MatAXPY(MatDeltaFwdSensip,-1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
773:     }
774:   } else { /* 1-stage method */
775:     th->shift = 0.0;
776:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
777:     MatMatMult(J,ts->mat_sensip,MAT_REUSE_MATRIX,PETSC_DEFAULT,&MatDeltaFwdSensip);
778:     MatScale(MatDeltaFwdSensip,-1.);

780:     /* Add the f_p forcing terms */
781:     if (ts->Jacp) {
782:       TSComputeIJacobianP(ts,th->stage_time,th->X,th->Xdot,th->shift,ts->Jacp,PETSC_FALSE);
783:       MatAXPY(MatDeltaFwdSensip,-1.,ts->Jacp,SUBSET_NONZERO_PATTERN);
784:     }
785:   }

787:   /* Build LHS */
788:   th->shift = 1/(th->Theta*th->time_step);
789:   if (th->endpoint) {
790:     TSComputeIJacobian(ts,th->stage_time,ts->vec_sol,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
791:   } else {
792:     TSComputeIJacobian(ts,th->stage_time,th->X,th->Xdot,th->shift,J,Jpre,PETSC_FALSE);
793:   }
794:   KSPSetOperators(ksp,J,Jpre);

796:   /*
797:     Evaluate the first stage of integral gradients with the 2-stage method:
798:     drdu|t_n*S(t_n) + drdp|t_n
799:     This is done before the linear solve because the sensitivity variable S(t_n) will be propagated to S(t_{n+1})
800:   */
801:   if (th->endpoint) { /* 2-stage method only */
802:     if (quadts && quadts->mat_sensip) {
803:       TSComputeRHSJacobian(quadts,th->ptime,th->X0,quadJ,NULL);
804:       TSComputeRHSJacobianP(quadts,th->ptime,th->X0,quadJp);
805:       MatTransposeMatMult(ts->mat_sensip,quadJ,MAT_REUSE_MATRIX,PETSC_DEFAULT,&th->MatIntegralSensipTemp);
806:       MatAXPY(th->MatIntegralSensipTemp,1,quadJp,SAME_NONZERO_PATTERN);
807:       MatAXPY(quadts->mat_sensip,th->time_step*(1.-th->Theta),th->MatIntegralSensipTemp,SAME_NONZERO_PATTERN);
808:     }
809:   }

811:   /* Solve the tangent linear equation for forward sensitivities to parameters */
812:   for (ntlm=0; ntlm<th->num_tlm; ntlm++) {
813:     KSPConvergedReason kspreason;
814:     MatDenseGetColumn(MatDeltaFwdSensip,ntlm,&barr);
815:     VecPlaceArray(VecDeltaFwdSensipCol,barr);
816:     if (th->endpoint) {
817:       MatDenseGetColumn(ts->mat_sensip,ntlm,&xarr);
818:       VecPlaceArray(ts->vec_sensip_col,xarr);
819:       KSPSolve(ksp,VecDeltaFwdSensipCol,ts->vec_sensip_col);
820:       VecResetArray(ts->vec_sensip_col);
821:       MatDenseRestoreColumn(ts->mat_sensip,&xarr);
822:     } else {
823:       KSPSolve(ksp,VecDeltaFwdSensipCol,VecDeltaFwdSensipCol);
824:     }
825:     KSPGetConvergedReason(ksp,&kspreason);
826:     if (kspreason < 0) {
827:       ts->reason = TSFORWARD_DIVERGED_LINEAR_SOLVE;
828:       PetscInfo2(ts,"Step=%D, %Dth tangent linear solve, linear solve fails, stopping tangent linear solve\n",ts->steps,ntlm);
829:     }
830:     VecResetArray(VecDeltaFwdSensipCol);
831:     MatDenseRestoreColumn(MatDeltaFwdSensip,&barr);
832:   }


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_step,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_step*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;

868:   if (ns) *ns = 1;
869:   if (stagesensip) *stagesensip = th->endpoint ? &(th->MatFwdSensip0) : &(th->MatDeltaFwdSensip);
870:   return(0);
871: }

873: /*------------------------------------------------------------*/
874: static PetscErrorCode TSReset_Theta(TS ts)
875: {
876:   TS_Theta       *th = (TS_Theta*)ts->data;

880:   VecDestroy(&th->X);
881:   VecDestroy(&th->Xdot);
882:   VecDestroy(&th->X0);
883:   VecDestroy(&th->affine);

885:   VecDestroy(&th->vec_sol_prev);
886:   VecDestroy(&th->vec_lte_work);

888:   VecDestroy(&th->VecCostIntegral0);
889:   return(0);
890: }

892: static PetscErrorCode TSAdjointReset_Theta(TS ts)
893: {
894:   TS_Theta       *th = (TS_Theta*)ts->data;

898:   VecDestroyVecs(ts->numcost,&th->VecsDeltaLam);
899:   VecDestroyVecs(ts->numcost,&th->VecsDeltaMu);
900:   VecDestroyVecs(ts->numcost,&th->VecsDeltaLam2);
901:   VecDestroyVecs(ts->numcost,&th->VecsDeltaMu2);
902:   VecDestroyVecs(ts->numcost,&th->VecsSensiTemp);
903:   VecDestroyVecs(ts->numcost,&th->VecsSensi2Temp);
904:   return(0);
905: }

907: static PetscErrorCode TSDestroy_Theta(TS ts)
908: {

912:   TSReset_Theta(ts);
913:   if (ts->dm) {
914:     DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
915:     DMSubDomainHookRemove(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);
916:   }
917:   PetscFree(ts->data);
918:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",NULL);
919:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",NULL);
920:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",NULL);
921:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",NULL);
922:   return(0);
923: }

925: /*
926:   This defines the nonlinear equation that is to be solved with SNES
927:   G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0
928: */
929: static PetscErrorCode SNESTSFormFunction_Theta(SNES snes,Vec x,Vec y,TS ts)
930: {
931:   TS_Theta       *th = (TS_Theta*)ts->data;
933:   Vec            X0,Xdot;
934:   DM             dm,dmsave;
935:   PetscReal      shift = th->shift;

938:   SNESGetDM(snes,&dm);
939:   /* When using the endpoint variant, this is actually 1/Theta * Xdot */
940:   TSThetaGetX0AndXdot(ts,dm,&X0,&Xdot);
941:   if (x != X0) {
942:     VecAXPBYPCZ(Xdot,-shift,shift,0,X0,x);
943:   } else {
944:     VecZeroEntries(Xdot);
945:   }
946:   /* DM monkey-business allows user code to call TSGetDM() inside of functions evaluated on levels of FAS */
947:   dmsave = ts->dm;
948:   ts->dm = dm;
949:   TSComputeIFunction(ts,th->stage_time,x,Xdot,y,PETSC_FALSE);
950:   ts->dm = dmsave;
951:   TSThetaRestoreX0AndXdot(ts,dm,&X0,&Xdot);
952:   return(0);
953: }

955: static PetscErrorCode SNESTSFormJacobian_Theta(SNES snes,Vec x,Mat A,Mat B,TS ts)
956: {
957:   TS_Theta       *th = (TS_Theta*)ts->data;
959:   Vec            Xdot;
960:   DM             dm,dmsave;
961:   PetscReal      shift = th->shift;

964:   SNESGetDM(snes,&dm);
965:   /* Xdot has already been computed in SNESTSFormFunction_Theta (SNES guarantees this) */
966:   TSThetaGetX0AndXdot(ts,dm,NULL,&Xdot);

968:   dmsave = ts->dm;
969:   ts->dm = dm;
970:   TSComputeIJacobian(ts,th->stage_time,x,Xdot,shift,A,B,PETSC_FALSE);
971:   ts->dm = dmsave;
972:   TSThetaRestoreX0AndXdot(ts,dm,NULL,&Xdot);
973:   return(0);
974: }

976: static PetscErrorCode TSForwardSetUp_Theta(TS ts)
977: {
978:   TS_Theta       *th = (TS_Theta*)ts->data;
979:   TS             quadts = ts->quadraturets;

983:   /* combine sensitivities to parameters and sensitivities to initial values into one array */
984:   th->num_tlm = ts->num_parameters;
985:   MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatDeltaFwdSensip);
986:   if (quadts && quadts->mat_sensip) {
987:     MatDuplicate(quadts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatIntegralSensipTemp);
988:     MatDuplicate(quadts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatIntegralSensip0);
989:   }
990:   /* backup sensitivity results for roll-backs */
991:   MatDuplicate(ts->mat_sensip,MAT_DO_NOT_COPY_VALUES,&th->MatFwdSensip0);

993:   VecDuplicate(ts->vec_sol,&th->VecDeltaFwdSensipCol);
994:   return(0);
995: }

997: static PetscErrorCode TSForwardReset_Theta(TS ts)
998: {
999:   TS_Theta       *th = (TS_Theta*)ts->data;
1000:   TS             quadts = ts->quadraturets;

1004:   if (quadts && quadts->mat_sensip) {
1005:     MatDestroy(&th->MatIntegralSensipTemp);
1006:     MatDestroy(&th->MatIntegralSensip0);
1007:   }
1008:   VecDestroy(&th->VecDeltaFwdSensipCol);
1009:   MatDestroy(&th->MatDeltaFwdSensip);
1010:   MatDestroy(&th->MatFwdSensip0);
1011:   return(0);
1012: }

1014: static PetscErrorCode TSSetUp_Theta(TS ts)
1015: {
1016:   TS_Theta       *th = (TS_Theta*)ts->data;
1017:   TS             quadts = ts->quadraturets;
1018:   PetscBool      match;

1022:   if (!th->VecCostIntegral0 && quadts && ts->costintegralfwd) { /* back up cost integral */
1023:     VecDuplicate(quadts->vec_sol,&th->VecCostIntegral0);
1024:   }
1025:   if (!th->X) {
1026:     VecDuplicate(ts->vec_sol,&th->X);
1027:   }
1028:   if (!th->Xdot) {
1029:     VecDuplicate(ts->vec_sol,&th->Xdot);
1030:   }
1031:   if (!th->X0) {
1032:     VecDuplicate(ts->vec_sol,&th->X0);
1033:   }
1034:   if (th->endpoint) {
1035:     VecDuplicate(ts->vec_sol,&th->affine);
1036:   }

1038:   th->order = (th->Theta == 0.5) ? 2 : 1;
1039:   th->shift = 1/(th->Theta*ts->time_step);

1041:   TSGetDM(ts,&ts->dm);
1042:   DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSTheta,DMRestrictHook_TSTheta,ts);
1043:   DMSubDomainHookAdd(ts->dm,DMSubDomainHook_TSTheta,DMSubDomainRestrictHook_TSTheta,ts);

1045:   TSGetAdapt(ts,&ts->adapt);
1046:   TSAdaptCandidatesClear(ts->adapt);
1047:   PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&match);
1048:   if (!match) {
1049:     VecDuplicate(ts->vec_sol,&th->vec_sol_prev);
1050:     VecDuplicate(ts->vec_sol,&th->vec_lte_work);
1051:   }
1052:   TSGetSNES(ts,&ts->snes);
1053:   return(0);
1054: }

1056: /*------------------------------------------------------------*/

1058: static PetscErrorCode TSAdjointSetUp_Theta(TS ts)
1059: {
1060:   TS_Theta       *th = (TS_Theta*)ts->data;

1064:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam);
1065:   VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsSensiTemp);
1066:   if (ts->vecs_sensip) {
1067:     VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&th->VecsDeltaMu);
1068:   }
1069:   if (ts->vecs_sensi2) {
1070:     VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&th->VecsDeltaLam2);
1071:     VecDuplicateVecs(ts->vecs_sensi2[0],ts->numcost,&th->VecsSensi2Temp);
1072:     /* hack ts to make implicit TS solver work when users provide only explicit versions of callbacks (RHSFunction,RHSJacobian,RHSHessian etc.) */
1073:     if (!ts->ihessianproduct_fuu) ts->vecs_fuu = ts->vecs_guu;
1074:     if (!ts->ihessianproduct_fup) ts->vecs_fup = ts->vecs_gup;
1075:   }
1076:   if (ts->vecs_sensi2p) {
1077:     VecDuplicateVecs(ts->vecs_sensi2p[0],ts->numcost,&th->VecsDeltaMu2);
1078:     /* hack ts to make implicit TS solver work when users provide only explicit versions of callbacks (RHSFunction,RHSJacobian,RHSHessian etc.) */
1079:     if (!ts->ihessianproduct_fpu) ts->vecs_fpu = ts->vecs_gpu;
1080:     if (!ts->ihessianproduct_fpp) ts->vecs_fpp = ts->vecs_gpp;
1081:   }
1082:   return(0);
1083: }

1085: static PetscErrorCode TSSetFromOptions_Theta(PetscOptionItems *PetscOptionsObject,TS ts)
1086: {
1087:   TS_Theta       *th = (TS_Theta*)ts->data;

1091:   PetscOptionsHead(PetscOptionsObject,"Theta ODE solver options");
1092:   {
1093:     PetscOptionsReal("-ts_theta_theta","Location of stage (0<Theta<=1)","TSThetaSetTheta",th->Theta,&th->Theta,NULL);
1094:     PetscOptionsBool("-ts_theta_endpoint","Use the endpoint instead of midpoint form of the Theta method","TSThetaSetEndpoint",th->endpoint,&th->endpoint,NULL);
1095:     PetscOptionsBool("-ts_theta_initial_guess_extrapolate","Extrapolate stage initial guess from previous solution (sometimes unstable)","TSThetaSetExtrapolate",th->extrapolate,&th->extrapolate,NULL);
1096:   }
1097:   PetscOptionsTail();
1098:   return(0);
1099: }

1101: static PetscErrorCode TSView_Theta(TS ts,PetscViewer viewer)
1102: {
1103:   TS_Theta       *th = (TS_Theta*)ts->data;
1104:   PetscBool      iascii;

1108:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1109:   if (iascii) {
1110:     PetscViewerASCIIPrintf(viewer,"  Theta=%g\n",(double)th->Theta);
1111:     PetscViewerASCIIPrintf(viewer,"  Extrapolation=%s\n",th->extrapolate ? "yes" : "no");
1112:   }
1113:   return(0);
1114: }

1116: static PetscErrorCode TSThetaGetTheta_Theta(TS ts,PetscReal *theta)
1117: {
1118:   TS_Theta *th = (TS_Theta*)ts->data;

1121:   *theta = th->Theta;
1122:   return(0);
1123: }

1125: static PetscErrorCode TSThetaSetTheta_Theta(TS ts,PetscReal theta)
1126: {
1127:   TS_Theta *th = (TS_Theta*)ts->data;

1130:   if (theta <= 0 || 1 < theta) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Theta %g not in range (0,1]",(double)theta);
1131:   th->Theta = theta;
1132:   th->order = (th->Theta == 0.5) ? 2 : 1;
1133:   return(0);
1134: }

1136: static PetscErrorCode TSThetaGetEndpoint_Theta(TS ts,PetscBool *endpoint)
1137: {
1138:   TS_Theta *th = (TS_Theta*)ts->data;

1141:   *endpoint = th->endpoint;
1142:   return(0);
1143: }

1145: static PetscErrorCode TSThetaSetEndpoint_Theta(TS ts,PetscBool flg)
1146: {
1147:   TS_Theta *th = (TS_Theta*)ts->data;

1150:   th->endpoint = flg;
1151:   return(0);
1152: }

1154: #if defined(PETSC_HAVE_COMPLEX)
1155: static PetscErrorCode TSComputeLinearStability_Theta(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi)
1156: {
1157:   PetscComplex   z   = xr + xi*PETSC_i,f;
1158:   TS_Theta       *th = (TS_Theta*)ts->data;
1159:   const PetscReal one = 1.0;

1162:   f   = (one + (one - th->Theta)*z)/(one - th->Theta*z);
1163:   *yr = PetscRealPartComplex(f);
1164:   *yi = PetscImaginaryPartComplex(f);
1165:   return(0);
1166: }
1167: #endif

1169: static PetscErrorCode TSGetStages_Theta(TS ts,PetscInt *ns,Vec **Y)
1170: {
1171:   TS_Theta     *th = (TS_Theta*)ts->data;

1174:   if (ns) *ns = 1;
1175:   if (Y)  *Y  = th->endpoint ? &(th->X0) : &(th->X);
1176:   return(0);
1177: }

1179: /* ------------------------------------------------------------ */
1180: /*MC
1181:       TSTHETA - DAE solver using the implicit Theta method

1183:    Level: beginner

1185:    Options Database:
1186: +  -ts_theta_theta <Theta> - Location of stage (0<Theta<=1)
1187: .  -ts_theta_endpoint <flag> - Use the endpoint (like Crank-Nicholson) instead of midpoint form of the Theta method
1188: -  -ts_theta_initial_guess_extrapolate <flg> - Extrapolate stage initial guess from previous solution (sometimes unstable)

1190:    Notes:
1191: $  -ts_type theta -ts_theta_theta 1.0 corresponds to backward Euler (TSBEULER)
1192: $  -ts_type theta -ts_theta_theta 0.5 corresponds to the implicit midpoint rule
1193: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint corresponds to Crank-Nicholson (TSCN)

1195:    This method can be applied to DAE.

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

1199: .vb
1200:   Theta | Theta
1201:   -------------
1202:         |  1
1203: .ve

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

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

1209: .vb
1210:   0 | 0         0
1211:   1 | 1-Theta   Theta
1212:   -------------------
1213:     | 1-Theta   Theta
1214: .ve

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

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

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

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

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

1226: M*/
1227: PETSC_EXTERN PetscErrorCode TSCreate_Theta(TS ts)
1228: {
1229:   TS_Theta       *th;

1233:   ts->ops->reset           = TSReset_Theta;
1234:   ts->ops->adjointreset    = TSAdjointReset_Theta;
1235:   ts->ops->destroy         = TSDestroy_Theta;
1236:   ts->ops->view            = TSView_Theta;
1237:   ts->ops->setup           = TSSetUp_Theta;
1238:   ts->ops->adjointsetup    = TSAdjointSetUp_Theta;
1239:   ts->ops->adjointreset    = TSAdjointReset_Theta;
1240:   ts->ops->step            = TSStep_Theta;
1241:   ts->ops->interpolate     = TSInterpolate_Theta;
1242:   ts->ops->evaluatewlte    = TSEvaluateWLTE_Theta;
1243:   ts->ops->rollback        = TSRollBack_Theta;
1244:   ts->ops->setfromoptions  = TSSetFromOptions_Theta;
1245:   ts->ops->snesfunction    = SNESTSFormFunction_Theta;
1246:   ts->ops->snesjacobian    = SNESTSFormJacobian_Theta;
1247: #if defined(PETSC_HAVE_COMPLEX)
1248:   ts->ops->linearstability = TSComputeLinearStability_Theta;
1249: #endif
1250:   ts->ops->getstages       = TSGetStages_Theta;
1251:   ts->ops->adjointstep     = TSAdjointStep_Theta;
1252:   ts->ops->adjointintegral = TSAdjointCostIntegral_Theta;
1253:   ts->ops->forwardintegral = TSForwardCostIntegral_Theta;
1254:   ts->default_adapt_type   = TSADAPTNONE;

1256:   ts->ops->forwardsetup     = TSForwardSetUp_Theta;
1257:   ts->ops->forwardreset     = TSForwardReset_Theta;
1258:   ts->ops->forwardstep      = TSForwardStep_Theta;
1259:   ts->ops->forwardgetstages = TSForwardGetStages_Theta;

1261:   ts->usessnes = PETSC_TRUE;

1263:   PetscNewLog(ts,&th);
1264:   ts->data = (void*)th;

1266:   th->VecsDeltaLam    = NULL;
1267:   th->VecsDeltaMu     = NULL;
1268:   th->VecsSensiTemp   = NULL;
1269:   th->VecsSensi2Temp  = NULL;

1271:   th->extrapolate = PETSC_FALSE;
1272:   th->Theta       = 0.5;
1273:   th->order       = 2;
1274:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetTheta_C",TSThetaGetTheta_Theta);
1275:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetTheta_C",TSThetaSetTheta_Theta);
1276:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaGetEndpoint_C",TSThetaGetEndpoint_Theta);
1277:   PetscObjectComposeFunction((PetscObject)ts,"TSThetaSetEndpoint_C",TSThetaSetEndpoint_Theta);
1278:   return(0);
1279: }

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

1284:   Not Collective

1286:   Input Parameter:
1287: .  ts - timestepping context

1289:   Output Parameter:
1290: .  theta - stage abscissa

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

1295:   Level: Advanced

1297: .seealso: TSThetaSetTheta()
1298: @*/
1299: PetscErrorCode  TSThetaGetTheta(TS ts,PetscReal *theta)
1300: {

1306:   PetscUseMethod(ts,"TSThetaGetTheta_C",(TS,PetscReal*),(ts,theta));
1307:   return(0);
1308: }

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

1313:   Not Collective

1315:   Input Parameter:
1316: +  ts - timestepping context
1317: -  theta - stage abscissa

1319:   Options Database:
1320: .  -ts_theta_theta <theta>

1322:   Level: Intermediate

1324: .seealso: TSThetaGetTheta()
1325: @*/
1326: PetscErrorCode  TSThetaSetTheta(TS ts,PetscReal theta)
1327: {

1332:   PetscTryMethod(ts,"TSThetaSetTheta_C",(TS,PetscReal),(ts,theta));
1333:   return(0);
1334: }

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

1339:   Not Collective

1341:   Input Parameter:
1342: .  ts - timestepping context

1344:   Output Parameter:
1345: .  endpoint - PETSC_TRUE when using the endpoint variant

1347:   Level: Advanced

1349: .seealso: TSThetaSetEndpoint(), TSTHETA, TSCN
1350: @*/
1351: PetscErrorCode TSThetaGetEndpoint(TS ts,PetscBool *endpoint)
1352: {

1358:   PetscUseMethod(ts,"TSThetaGetEndpoint_C",(TS,PetscBool*),(ts,endpoint));
1359:   return(0);
1360: }

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

1365:   Not Collective

1367:   Input Parameter:
1368: +  ts - timestepping context
1369: -  flg - PETSC_TRUE to use the endpoint variant

1371:   Options Database:
1372: .  -ts_theta_endpoint <flg>

1374:   Level: Intermediate

1376: .seealso: TSTHETA, TSCN
1377: @*/
1378: PetscErrorCode TSThetaSetEndpoint(TS ts,PetscBool flg)
1379: {

1384:   PetscTryMethod(ts,"TSThetaSetEndpoint_C",(TS,PetscBool),(ts,flg));
1385:   return(0);
1386: }

1388: /*
1389:  * TSBEULER and TSCN are straightforward specializations of TSTHETA.
1390:  * The creation functions for these specializations are below.
1391:  */

1393: static PetscErrorCode TSSetUp_BEuler(TS ts)
1394: {
1395:   TS_Theta       *th = (TS_Theta*)ts->data;

1399:   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");
1400:   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");
1401:   TSSetUp_Theta(ts);
1402:   return(0);
1403: }

1405: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
1406: {
1408:   return(0);
1409: }

1411: /*MC
1412:       TSBEULER - ODE solver using the implicit backward Euler method

1414:   Level: beginner

1416:   Notes:
1417:   TSBEULER is equivalent to TSTHETA with Theta=1.0

1419: $  -ts_type theta -ts_theta_theta 1.0

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

1423: M*/
1424: PETSC_EXTERN PetscErrorCode TSCreate_BEuler(TS ts)
1425: {

1429:   TSCreate_Theta(ts);
1430:   TSThetaSetTheta(ts,1.0);
1431:   TSThetaSetEndpoint(ts,PETSC_FALSE);
1432:   ts->ops->setup = TSSetUp_BEuler;
1433:   ts->ops->view  = TSView_BEuler;
1434:   return(0);
1435: }

1437: static PetscErrorCode TSSetUp_CN(TS ts)
1438: {
1439:   TS_Theta       *th = (TS_Theta*)ts->data;

1443:   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");
1444:   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");
1445:   TSSetUp_Theta(ts);
1446:   return(0);
1447: }

1449: static PetscErrorCode TSView_CN(TS ts,PetscViewer viewer)
1450: {
1452:   return(0);
1453: }

1455: /*MC
1456:       TSCN - ODE solver using the implicit Crank-Nicolson method.

1458:   Level: beginner

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

1463: $  -ts_type theta -ts_theta_theta 0.5 -ts_theta_endpoint

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

1467: M*/
1468: PETSC_EXTERN PetscErrorCode TSCreate_CN(TS ts)
1469: {

1473:   TSCreate_Theta(ts);
1474:   TSThetaSetTheta(ts,0.5);
1475:   TSThetaSetEndpoint(ts,PETSC_TRUE);
1476:   ts->ops->setup = TSSetUp_CN;
1477:   ts->ops->view  = TSView_CN;
1478:   return(0);
1479: }